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Chapter  1 

Introduction 


The  main  objective  of  the  work  described  in  this  report  is  to  better  understand  the  flow  physics 
of  aircraft  wings  undergoing  highly  unsteady  maneuvers.  Reduced-order  models  play  a  central  role 
in  this  study,  both  to  elucidate  the  overall  dynamical  mechanisms  behind  various  flow  phenomena 
(such  as  dynamic  stall  and  vortex  shedding),  and  ultimately  to  guide  flight  control  design  for 
vehicles  for  which  these  unsteady  phenomena  are  important. 

Unsteady  phenomena  are  of  increasing  interest  to  the  Air  Force,  as  lightweight  unmanned  air 
vehicles  become  more  prevalent.  With  increasingly  smaller  and  lighter  vehicles  envisioned  in  the 
future,  understanding  unsteady  aerodynamics  will  become  even  more  important,  in  order  to  design 
control  systems  that  can  respond  to  severe  gusts,  or  perform  highly  agile  maneuvers. 

Our  approach  builds  upon  recent  advances  in  understanding  the  dynamics  of  these  unsteady 
flows,  and  uses  state-of-the-art  techniques,  both  for  measuring  these  phenomena  in  experiments 
(using  an  unsteady  wind  tunnel  at  IIT),  and  for  analyzing  the  data  and  developing  reduced-order 
models  (using  techniques  such  as  the  Eigensystem  Realization  Algorithm  and  variants  of  Dynamic 
Mode  Decomposition,  being  developed  at  Princeton). 

The  studying  of  aerodynamics  and  fluid  mechanics  predate  the  birth  of  the  scientific  method 
itself.  From  the  initial  musings  of  Aristotle  and  Archimedes,  to  the  sketches  of  Leonardo  Da  Vinci, 
to  George  Caley’s  conception  of  a  modern  airplane  configuration  in  1799.  Yet,  in  spite  of  this  long 
history,  and  the  fact  that  the  underlying  governing  equations  of  viscous  fluid  flow  being  known  for 
hundreds  of  years,  there  is  still  much  to  learn.  For  example,  in  many  respects  we  are  still  catching 
up  to  the  mastery  exhibited  by  biological  swimmers  and  fliers,  in  terms  of  maneuverability  and 
efficiency.  This  “unfinished  business”  is  not  by  any  means  due  to  a  lack  of  imperative.  Indeed, 
reductions  in  aerodynamic  drag  on  cars,  trucks,  airplanes  and  ships  can  lead  to  billions  of  dollars  in 
fuel  savings,  and  significant  reductions  of  CO2  emissions,  to  speak  nothing  of  the  enhanced  safety 
and  maneuverability  that  would  come  with,  for  example,  a  more  comprehensive  understanding  of 
the  dynamics  of  aircraft  in  deep  stall. 

With  advances  in  both  experimental  techniques  and  equipment,  and  computational  power  and 
storage  capacity,  researchers  in  fluid  dynamics  can  now  generate  more  high-fidelity  data  than  ever 
before.  This  is  crucial  to  advancing  the  field  as  a  whole,  since  all  but  the  most  idealized,  simple 
systems  cannot  be  fully  understood  by  analytical  deductions  alone.  Moreover,  advances  in  mate¬ 
rials,  manufacturing,  and  the  miniaturization  of  processors  allow  for  new  questions  to  become  of 
practical  interest.  This  is  not  to  suggest  that  the  field  of  fluid  dynamics  and  aerodynamics  are  now 
confined  to  the  realms  of  data  science.  Indeed,  progress  is  made  by  combining  the  insight  attained 
from  study  of  the  fundamental  equations  with  the  additional  dimension  of  large  data.  Indeed,  ob¬ 
taining  the  right  data  relies  on  understanding  of  the  physical  system.  Beyond  this,  however,  data 
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collection  should  also  be  informed  by  a  proper  understanding  of  the  capabilities  and  limitations  of 
the  algorithmic  tools  that  are  required  to  assist  in  data  analysis. 

The  presence  of  increasingly  large  data  sets  necessitates  the  use  of  such  post-processing  tech¬ 
niques  that  are  able  to  extract  tractable  and  physically  relevant  information  from  the  data.  This 
report  considers  methods  to  extract  pertinent  information  and  models  from  data  obtained  from  flu¬ 
ids  simulations  and  experiments.  While  the  focus  is  on  unsteady  aerodynamic  systems,  specifically 
pitching  airfoils,  the  techniques  used  and  developed  are  applicable  to  a  wider  range  of  applications, 
both  within  and  outside  of  fluids  systems. 

1.1  Unsteady  aerodynamics 

The  flight  of  small,  highly  maneuverable  aircraft,  whether  biological  or  manmade,  is  greatly  im¬ 
pacted  by  unsteady  aerodynamic  effects,  which  can  be  either  beneficial  or  detrimental  to  flight. 
Accurate  understanding  of  such  effects  can  allow  for  the  design  of  aircraft  that  are  more  efficient, 
responsive,  and  robust. 

The  need  to  account  for  unsteady  effects  has  been  recognized  since  soon  after  the  breakthrough  of 
powered  manmade  flight,  in  the  classical  works  of  Wagner  [155],  Theodorsen  [142],  and  Garrick  [50]. 
Indeed,  many  failed  attempts  at  flight  can  probably  be  attributed  to  a  severe  lack  of  understanding 
of  how  to  utilize  such  effects.  These  classical  models  give  significant  insight  into  the  fundamental 
flow  physics  associated  with  unsteady  flight,  such  as  relative  contributions  to  lift  of  the  added  mass, 
quasi-steady  bound  circulation,  and  wake  vortices.  For  example,  the  Theodorsen  model  gives  the 
relationship  between  the  airfoil  kinematics  (a,  a,  and  a  and  the  lift  coefficient  by: 

Cl  =  |  (a  -  +  2vr  Q  -  C(k),  (1-1.1) 

where  a  is  a  parameter  that  defines  the  pitch  axes,  with  a  =  —  1  and  +1  corresponding  to  pitching 
about  the  leading  and  trailing  edge  of  the  airfoil,  respectively,  k  =  ^  is  the  reduced  frequency, 
and  C(k)  is  the  Theodorsen  function,  which  governs  the  interaction  between  the  shed  vorticity  in 
the  wake  and  the  circulatory  lift  force.  While  such  models  can  be  quantitatively  accurate  for  cases 
of  attached  flow  where  viscous  effects  are  negligible,  they  quickly  lose  validity  when  dealing  with 
separated  flows,  which  are  often  encountered  in  the  extreme  motions  that  are  possible  for  birds, 
insects,  and  micro  and  unmanned  aerial  vehicles  (MAV  and  UAV).  It  is  precisely  in  these  extreme 
cases  that  accurate  predictive  models  are  essential  to  prevent  catastrophic  failure  and  ensure  ongoing 
successful  flight.  While  more  accurate  predictions  can  be  attained  from  high-fidelity  simulations, 
the  computational  cost  typically  prohibits  the  direct  use  of  such  simulations  for  real-time  prediction 
and  control. 

Biological  examples  such  as  insects  [18,  120,  156]  and  birds  [153]  have  seemingly  evolved  to  take 
advantage  of  the  high  transient  lift  force  that  can  be  generated  due  to  the  formation  of  a  leading 
edge  vortex  (LEV)  during  rapid  pitch-  up  motion,  for  example.  While  these  give  motivating  ex¬ 
amples  of  the  advantages  of  accurate  understanding  of  unsteady  aerodynamic  effects,  the  preferred 
wing  kinematics  arising  from  evolution  is  highly  specific  and  coupled  to  the  geometry  and  other 
physiological  features  of  the  animal.  Indeed,  the  characteristics  of  unsteady  aerodynamic  effects, 
particularly  for  separated  flows,  seem  to  be  quite  sensitive  to  both  the  geometry  [79]  and  Reynolds 
number  [166]  of  the  airfoils.  Studies  into  low  Reynolds  number  flow  over  stationary  [166,  4,  165] 
and  pitching  [154,  1,  131,  30,  80]  symmetric  airfoils  have  revealed,  for  example,  complex  effects 
associated  with  the  stability  and  separation  of  the  suction  surface  boundary  layer,  which  are  again 
highly  sensitive  to  Reynolds  numbers.  These  observations  motivate  the  development  of  general 
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modeling  procedures  that  can  be  easily  applied  to  a  range  of  parameter  cases.  In  addition,  it  is 
desirable  for  such  methods  to  be  sufficiently  general  such  that  they  can  be  applied  to  more  realistic 
aircraft  configurations,  rather  than  just  airfoils.  As  an  example,  such  data  driven  modeling  was 
considered  for  the  case  of  accurate  prediction  and  control  of  lift  for  a  low  Reynolds  number  pitching 
airfoil  [22,  23],  using  the  eigensystem  realization  algorithm  [72]  (ERA)  and  observer /Kalman  filter 
identification  [74]  (OKID).  There  has  also  been  a  significant  amount  of  work  in  terms  of  nonlinear 
modeling,  ranging  from  low  order  state-  space  models  formulated  from  theoretical  considerations 
[54],  to  Volterra  series  models  that  have  been  used  to  model  a  range  of  unsteady  aerodynamic  and 
aeroelastic  phenomena  [132,  85,  11]. 


1.2  Data-driven  modeling 

With  advances  in  both  experimental  techniques  and  equipment,  and  computational  power  and 
storage  capacity,  researchers  in  fluid  dynamics  can  now  generate  more  high-fidelity  data  than  ever 
before.  The  presence  of  increasingly  large  data  sets  calls  for  appropriate  data  analysis  techniques, 
that  are  able  to  extract  tractable  and  physically  relevant  information  from  the  data.  In  particular,  a 
much-desired  goal  in  fluid  mechanics,  and  indeed  many  other  fields,  is  to  obtain  simple  models  that 
are  capable  of  predicting  the  behavior  of  seemingly  complex  systems.  Low- dimensional  models  can 
not  only  improve  our  fundamental  understanding  of  such  systems,  but  are  often  required  for  purpose 
of  efficient  and  accurate  prediction,  estimation  and  control.  Broadly  speaking,  one  can  obtain  low¬ 
dimensional  information  about  a  system  (whether  it  be  in  the  form  of  a  reduced-order  model,  or 
simply  spatial  modes  corresponding  to  certain  energetic  or  dynamic  characteristics)  in  numerous 
ways,  potentially  using  some  combination  of  data  collected  from  simulations  and  experiments,  and 
theoretical  knowledge  of  the  system,  such  as  the  governing  partial  differential  equations  (PDEs). 

Purely  data-driven  methods  can  include  those  developed  particularly  for  fluids  applications, 
such  as  the  dynamic  mode  decomposition  (DMD)  [125,  126],  or  those  which  are  appropriated  from 
other  communities,  such  as  the  eigensystem  realization  algorithm  (ERA)  [82,  72],  which  was  first 
applied  to  study  spacecraft  structures,  but  has  more  recently  been  appropriated  to  model  a  wide 
range  of  fluids  systems  [25,  2,  68,  69,  22,  23,  15,  67,  49]. 

Dynamic  mode  decomposition  allows  for  the  identification  and  analysis  of  dynamical  features  of 
time-evolving  fluid  flows,  using  data  obtained  from  either  experiments  or  simulations.  In  contrast 
to  other  data-driven  modal  decompositions  such  as  the  proper  orthogonal  decomposition  (POD), 
DMD  allows  for  spatial  modes  to  be  identified  that  can  be  directly  associated  with  characteristic 
frequencies  and  growth/decay  rates.  Following  its  conception,  DMD  was  quickly  shown  to  be 
useful  in  extracting  dynamical  features  in  both  experimental  and  numerical  data  [125,  126].  It 
has  subsequently  been  used  to  gain  dynamic  insight  on  a  wide  range  of  problems  arising  in  fluid 
mechanics  [e.g. ,  119,  128,  127,  100,  129,  43,  58,  70,  83,  93,  55,  123,  45]  and  other  fields  [e.g.,  59]. 

DMD  has  a  strong  connection  to  Koopman  operator  theory  [81,  88],  as  exposed  in  [119],  and 
further  reviewed  in  [89],  which  can  justify  its  use  in  analyzing  nonlinear  dynamical  systems.  Since 
its  original  formulation,  numerous  modifications  and  extensions  have  been  made  to  DMD.  Chen 
et  al.  [27]  highlights  the  connection  that  DMD  shares  with  traditional  Fourier  analysis,  as  well 
as  proposing  an  optimized  algorithm  that  recasts  DMD  as  an  optimal  dimensionality  reduction 
problem.  This  concept  of  finding  only  the  dynamically  important  modes  has  also  been  considered 
in  subsequent  works  of  [164]  and  [71].  All  of  these  works  are  motivated,  in  part,  by  the  fact  that 
by  default  DMD  will  output  as  many  modes  as  there  are  pairs  of  snapshots  (assuming  that  the 
length  of  the  snapshot  vector  is  greater  than  the  number  of  snapshots),  which  is  arbitrary  with 
respect  to  the  dynamical  system  under  consideration.  In  reality,  one  would  prefer  to  output  only 
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the  modes  and  eigenvalues  that  are  present  (or  dominant)  in  the  data.  When  the  data  is  corrupted 
by  noise  (as  will  always  be  the  case  to  some  degree,  especially  for  experimental  data),  this  process 
becomes  nontrivial,  since  noisy  data  might  have  a  numerical  rank  far  larger  than  the  dimension  of 
the  governing  dynamics  of  the  system.  Further  to  this,  one  cannot  expect  to  have  a  clean  partition 
into  modes  that  identify  true  dynamical  features,  and  those  which  consist  largely  of  noise. 

Simple  ways  of  achieving  this  objective  can  involve  either  first  projecting  the  data  onto  a  smaller 
dimensional  basis  (such  as  the  most  energetic  POD  modes)  before  applying  DMD,  or  by  choosing 
only  the  most  dynamically  important  DMD  modes  after  applying  DMD  to  the  full  data.  One 
can  also  truncate  the  data  to  a  dimension  larger  than  the  assumed  dimension  of  the  dynamics, 
and  then  apply  a  balanced  truncation  to  the  resulting  dynamical  system  to  obtain  the  desired 
reduced  order  model.  This  approach  is  sometimes  referred  to  as  over-specification  in  the  system 
identification  community  [see,  e.g.,  75].  Keeping  a  higher  dimension  of  data  than  that  of  the  assumed 
dynamics  can  be  particularly  important  for  input-output  systems  that  have  highly  energetic  modes 
that  are  not  strongly  observable  or  controllable  [112].  Ideally,  any  algorithm  that  restricts  the 
number  of  DMD  modes  that  are  computed  should  also  additionally  be  computationally  efficient. 
A  fast  method  to  perform  DMD  in  real  time  on  large  datasets  was  recently  proposed  in  [62], 
while  a  library  for  efficient  parallel  implementation  of  number  of  common  modal  decomposition 
and  system  identification  techniques  is  described  in  [16].  Sayadi  and  Schmid  [122]  also  gives  an 
explicit  implementation  of  DMD  for  parallelized  computation.  One  can  also  achieve  computational 
speedup  by  incorporating  efficient  methods  to  compute  singular  value  decompositions,  typically  the 
computational  bottleneck  in  DMD,  to  speed  up  the  computation  [] 

One  of  the  major  advantages  of  DMD  over  techniques  such  as  global  stability  analysis,  for  exam¬ 
ple,  is  that  it  can  be  applied  directly  to  data,  without  the  need  for  the  knowledge  or  construction  of 
the  system  matrix,  which  is  typically  not  available  for  experiments  [126].  For  this  reason,  analysis 
of  the  sensitivity  of  DMD  to  the  type  of  noise  typically  found  in  experimental  results  is  of  particular 
importance.  The  effects  of  noise  on  the  accuracy  of  the  DMD  procedure  was  systematically  inves¬ 
tigated  in  the  empirical  study  of  [44],  for  the  case  of  a  synthetic  waveform  inspired  by  canonical 
periodic  shear  flow  instabilities.  More  recently,  [101]  have  extended  this  type  of  analysis  to  more 
complex  data  with  multiple  frequencies,  as  might  be  found  in  typical  fluids  systems. 

An  notable  limitation  of  the  methods  mentioned  so  far  is  that  (when  considered  in  the  context 
of  data-driven  reduced-order  modeling  techniques)  they  are  linear,  in  the  sense  that  the  reduced 
order  model  that  is  identified  is  in  the  form  of  a  linear  system  of  ordinary  differential  equations 
(ODEs).  While  there  have  been  a  number  of  examples  of  nonlinear  data-driven  modeling  techniques 
used  in  fluids  applications  [85,  102,  104,  52,  11,  37,  77,  60,  39,  24],  their  widespread  use  has  been 
more  limited,  and  the  underlying  theory  is  less  established,  than  linear  techniques.  More  details 
concerning  the  application  of  data-driven  modeling  techniques  in  fluid  mechanics  can  be  found  in 
recent  review  articles  [21,  113]. 

Perhaps  the  most  common  method  to  obtain  a  nonlinear  reduced  order  model  for  fluids  systems 
comes  via  a  projection  of  the  governing  equations  onto  a  low-dimensional  basis  that  is  optimal  for 
capturing  the  energy  of  the  data,  i.e.,  the  proper  orthogonal  decomposition  (POD)  [86,  17,  66],  a 
procedure  referred  to  as  Galerkin  projection.  Galerkin  projection  (GP)  has  been  used  to  extract 
models  for  many  different  fluids  systems,  a  non-exhaustive  list  includes  flow  past  a  cylinder  at  low 
Reynolds  number  [41,  96,  95],  grooved  channels  [41]  the  wall  region  of  turbulent  boundary  layers 
[6,  103],  flat  plate  boundary  layers  [111],  turbulent  plane  Couette  flow  [91,  134],  turbulent  pipe  flow 
[19]  cavity  oscillations  [116,  115],  mixing  layers  [110,  150,  13],  and  compressible  flows  [118].  One 
significant  drawback  of  GP  models  is  that  they  ignore  modes  that  are  low  in  energy,  but  are  required 
for  the  dissipation  of  energy  in  the  full  system.  A  number  of  modifications  have  been  proposed  to 
address  this  concern,  as  well  as  other  issues  with  such  models.  Aubry  et  al.  [6]  and  Podvin  [103] 
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use  an  eddy  viscosity  term  that  accounts  for  energy  dissipated  into  unmodeled  modes,  Osth  et  al. 
[99]  investigate  a  hierachy  of  eddy  viscosity  formulations,  while  Wang  et  al.  [157,  158]  incorporate 
LES  closure  modeling  strategies.  Refs.  [35]  and  [36]  summarize  a  number  of  calibration  techniques 
that  can  be  used  to  improve  the  accuracy  of  Galerkin  models,  and  also  discuss  the  various  ways  in 
which  the  error  of  such  models  can  be  quantified.  Ref.  [12]  employs  a  subspace  rotation  technique 
to  stabilize  the  models,  which,  unlike  other  calibration  techniques,  maintains  consistency  with  the 
original  governing  equations.  Ref.  [13]  imposes  constraints  to  balance  the  turbulent  kinetic  energy 
of  the  resulting  model.  All  of  these  modifications  of  Galerkin  projection  increase  the  “data-driven” 
nature  of  the  method.  Ref.  [97]  gives  an  in-depth  summary  and  analysis  of  many  issues,  variations, 
progress,  and  open  problems  on  the  topic  of  Galerkin  projection  models,  while  [84]  gives  a  clear 
expository  introduction  of  the  main  ideas  in  Galerkin  projection,  with  examples. 

While  we  mentioned  above  that  DMD  could  be  classified  as  a  “linear”  method,  connections 
between  the  DMD  algorithm  and  the  Koopman  operator  [119,  89]  give  promise  that  it  can  ultimately 
be  used  to  model  and  understand  nonlinearities.  The  Koopman  operator  [81]  gives  a  means  of 
representing  a  finite-dimensional,  nonlinear  system  as  an  infinite  dimensional  linear  system,  and 
DMD  gives  a  finite-dimensional  approximation  to  this  operator. 

In  particular,  an  extension  of  DMD  that  potentially  allows  for  better  representation  of  non¬ 
linear  data  has  also  recently  been  proposed  [160],  and  although  the  computational  costs  increase 
dramatically  with  the  dimension  of  the  system,  a  kernel  method  described  in  [159]  reduces  the  cost 
to  be  comparable  to  standard  DMD. 

1.3  Organization  and  contributions 

Following  this  introductory  chapter,  Chapter  2  presents  a  summary  of  key  concepts  and  techniques 
in  data-driven  modeling  of  fluid  systems,  as  well  as  presenting  some  results  that,  besides  being  of 
some  independent  interest,  will  motivate  the  research  directions  taken  in  the  subsequent  chapters. 
Broadly  speaking,  Chapters  4  and  3  focus  specifically  on  the  application  pitching  airfoils,  while  5 
and  6  focus  on  data-driven  modeling  techniques,  particularly  extensions  and  improvements  to  DMD 
and  their  application  to  fluids  systems.  More  precisely,  the  report  includes  the  following  sections: 

Chapter  3  obtains  models  to  predict  the  pressures  and  forces  on  a  rapidly  picking  airfoil.  This  is 
one  of  the  first  applications  to  experimental  data  of  a  recently  developed  variant  of  DMD  to  allow  for 
the  identification  of  systems  that  contain  inputs  (e.g.,  systems  that  are  being  controlled  externally  in 
some  manner)  [105].  We  show  that  this  modeling  approach  is  convenient  for  constructing  “switched 
models”,  whereby  one  can  predict  the  behavior  of  a  nonlinear  system  by  switching  between  a 
family  of  linear  models.  In  particular,  the  “DMD  with  inputs”  method  gives  models  for  which  the 
coordinates  of  the  models  remain  consistent  with  each  other,  which  eliminates  complications  and 
ill-conditioning  that  has  been  observed  when  using  alternate  methods  [38] .  This  modeling  approach 
allows  for  the  formulation  of  a  switched  model  that  remains  accurate  over  a  wide  range  of  angles 
of  attack,  ranging  from  attached  to  fully  separated  flow.  The  experiments  and  subsequent  analysis 
were  performed  with  Nicole  Schiavone,  an  undergraduate  working  in  Prof.  Rowley’s  lab  in  summer 
2014.  Material  in  this  chapter  is  based  on  the  conference  paper: 

•  Scott  T  M  Dawson,  Nicole  K  Schiavone,  Clarence  W  Rowley,  and  David  R  Williams.  A  data- 
driven  modeling  framework  for  predicting  forces  and  pressures  on  a  rapidly  pitching  airfoil. 
In  45th  AIAA  Fluid  Dynamics  Conference,  page  2767,  2015. 

Chapter  4  explores  a  phenomena  that  is  identified  in  Chapter  2:  Namely  that  airfoils  under¬ 
going  low-amplitude  sinusoidal  pitching  motion  generate  enhanced  lift  when  pitching  occurs  at  a 
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preferred  frequencies.  A  systematic  parameter  sweep  over  the  pitching  frequency,  amplitude,  and 
base  angle  of  attack  is  conducted,  with  the  mean  and  frequency  content  of  the  forces  analyzed.  In 
addition,  the  flow  fields  are  studied  using  DMD.  Material  in  this  chapter  is  based  on  the  conference 
paper: 

•  Scott  TM  Dawson,  Daniel  C  Floryan,  Clarence  W  Rowley,  and  Maziar  S  Hemati.  Lift  en¬ 
hancement  of  high  angle  of  attack  airfoils  using  periodic  pitching.  In  54th  AIAA  Aerospace 
Sciences  Meeting,  page  2069,  2016. 

Chapter  5  shows  how  a  recently  developed  extension  to  DMD  can  be  utilized  to  obtain  nonlin¬ 
ear  reduced  order  models  for  fluids  systems.  We  modify  the  extended  DMD  algorithm  to  include  a 
Tikhonov  regularization  step,  which  is  found  to  give  improved  results  for  the  purposes  of  nonlinear 
system  identification.  The  method  is  demonstrated  on  the  canonical  example  of  flow  past  a  circular 
cylinder,  for  data  starting  near  the  unstable  equilibrium  solution  and  converging  to  the  periodic 
vortex  shedding  limit  cycle.  It  is  demonstrated  that  this  approach  can  be  superior  to  classical 
POD-Galerkin  projection,  particularly  in  cases  where  the  data  is  noisy,  is  from  a  limited  spatial 
region,  is  not  spatially  resolved,  or  is  only  collected  over  a  short  time  window.  Material  in  this 
chapter  is  contained  in  the  paper: 

•  Scott  T.  M.  Dawson  and  Clarence  W.  Rowley,  Nonlinear  reduced-order  models  of  fluids 
systems  using  extended  dynamic  mode  decomposition,  In  preparation  for  Theoretical  and 
Computational  Fluid  Dynamics ,  2016. 

Some  of  the  results  presented  in  Chapter  5  are  also  used  in  an  upcoming  review  paper: 

•  Clarence  W.  Rowley  and  Scott  T.  M.  Dawson.  Model  reduction  for  flow  analysis  and  control. 
Annual  Review  of  Fluid  Mechanics,  49(1),  2017. 

Chapter  6  analyzes  the  effect  of  noise  of  DMD.  As  well  as  giving  an  explanation  for  a  previously 
identifies  sensitivity  to  noisy  data,  three  variants  of  the  DMD  algorithm  are  proposed,  all  of  which 
perform  better  on  noisy  data.  This  work  was  performed  with  Maziar  Hemati  and  Matthew  Williams, 
and  is  published  in  the  following  paper: 

•  Scott  T.  M.  Dawson,  Maziar  S.  Hemati,  Matthew  O.  Williams,  and  Clarence  W.  Rowley. 
Characterizing  and  correcting  for  the  effect  of  sensor  noise  in  the  dynamic  mode  decomposi¬ 
tion.  Experiments  in  Fluids,  57(42):1?19,  2016. 

In  addition  to  the  publications  listed  above,  Chapter  2  uses  some  material  from  the  paper: 

•  Steven  L  Brunton,  Scott  T  M  Dawson,  and  Clarence  W  Rowley.  State-space  model  identifica¬ 
tion  and  feedback  control  of  unsteady  aerodynamic  forces.  Journal  of  Fluids  and  Structures, 
50:2537270,  2014. 

Effort  is  made  to  keep  notation  consistent  throughout  this  report;  however,  on  occasion  notation 
changes  between  chapters,  in  attempt  to  uphold  existing  conventions  in  the  relevant  fields.  While 
each  of  these  chapters  is  largely  self  contained,  we  present  in  Chapter  2  underlying  preliminaries 
that  are  broadly  relevant  across  all  sections  of  this  thesis. 

1.3.1  Publications  supported  by  this  grant 
Journal  papers: 
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•  Dawson  STM  and  Rowley  CW,  “Nonlinear  reduced-order  models  of  fluids  systems  using 
extended  dynamic  mode  decomposition,”  in  preparation. 

•  Hemati  MS,  Dawson  STM,  and  Rowley  CW,  “Parameter- Varying  Models  for  Unsteady  Aero¬ 
dynamic  Response  Prediction,”  submitted  to  AIAA  Journal. 

•  Rowley  CW  and  Dawson  STM,  “Model  reduction  for  flow  analysis  and  control,”  Annual 
Reviews  in  Fluid  Mechanics ,  to  appear,  Jan  2017. 

•  Monnier  B,  Williams  DR,  Weier  T,  and  Albrecht  T,  “Comparison  of  a  Separated  Flow  Re¬ 
sponse  to  Localized  and  Global  type  Disturbances,”  Experiments  in  Fluids ,  accepted  for 
publication,  Jun  2016. 

•  Williams  DR,  Reissner  F,  Greenblatt  D,  Muehler-Vahl  H,  and  Strangfeld  C,  “Modeling  Lift 
Hysteresis  on  Pitching  Airfoils  with  a  Modified  Goman-Khrabrov  Model,”  AIAA  Journal , 
accepted  for  publication,  Aug  2016. 

•  Dawson  STM,  Hemati  MS,  Williams  MO,  and  Rowley  CW,  “Characterizing  and  correcting 
for  the  effect  of  sensor  noise  in  the  dynamic  mode  decomposition,”  Experiments  in  Fluids, 
57:42,  2016. 

•  Williams  DR,  An  X,  Iliev  S,  King  R,  and  Reissner  F,  “Dynamic  hysteresis  control  of  lift  on 
a  pitching  wing,”  Experiments  in  Fluids,  56:1-12,  2015. 

•  An  X,  Grimaud  L,  and  Williams  D,  “Feedforward  Control  of  Lift  Hysteresis  during  Periodic 
and  Randomized  Pitching  Maneuvers,”  In:  King  R  (ed)  Active  Flow  Control  and  Combustion 
Control,  vol  127.  Springer,  Heidelberg,  pp  55-69. 

•  Tu  JH,  Rowley  CW,  Luchtenburg  DM,  Brunton  SL,  and  Kutz  JN,  “On  dynamic  mode  de¬ 
composition:  theory  and  applications,”  J.  Computational  Dynamics,  1(2):391-421,  Dec  2014. 

•  Hemati  MS,  Williams  MO,  and  Rowley  CW,  “Dynamic  mode  decomposition  for  large  and 
streaming  datasets,”  Physics  of  Fluids,  26(11):111781,  Nov  2014. 

•  Brunton  SL,  Dawson  STM,  and  Rowley  CW,  “State-space  model  and  feedback  control  of 
unsteady  aerodynamic  forces,”  Journal  of  Fluids  and  Structures,  50:253-270,  Oct  2014. 

•  Brunton  SL,  Rowley  CW,  and  Williams  DR,  “Reduced-order  unsteady  aerodynamic  models 
at  low  Reynolds  numbers,”  Journal  of  Fluid  Mechanics ,  724:203-233,  Sep  2013. 

Conference  papers: 

•  Dawson  STM,  Hemati  MS,  Floryan  D,  and  Rowley  CW,  “Lift  Enhancement  of  High  Angle  of 
Attack  Airfoils  using  Periodic  Pitching,”  54th  AIAA  Aerospace  Sciences  Meeting,  San  Diego, 
CA,  January  2015. 

•  Dawson  STM,  Schiavone  NK,  Rowley  CW,  and  Williams  DR,  “A  Data-Driven  Modeling 
Framework  for  Predicting  Forces  and  Pressures  on  a  Rapidly  Pitching  Airfoil,”  44th  AIAA 
Fluid  Mechanics  Conference,  Dallas,  TX,  June  2015. 

•  Hemati  MS,  Dawson  STM,  and  Rowley  CW,  “Unsteady  Aerodynamic  Response  Modeling: 
A  Parameter- Varying  Approach,”  53rd  AIAA  Aerospace  Sciences  Meeting,  Kissimmee,  FL, 
January  2014. 
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•  Monnier  B,  Williams  DR,  Weier  T,  and  Albrecht  T,  “Comparison  of  a  Separated  Flow  Re¬ 
sponse  to  Localized  and  Global-type  Disturbances,”  International  Journal  of  Flow  Control , 
AIAA  paper  2015-1056 

•  Williams  DR,  Reissner  F,  Greenblatt  D,  Muehler-Vahl  H,  and  Strangfeld  C,  “Modeling  Lift 
Hysteresis  with  a  Modified  Goman-Khrabrov  Model  on  Pitching  Airfoils,”  AIAA  Paper  2015- 
2631. 
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Chapter  2 

Background  and  motivating  results 


This  section  will  introduce  some  concepts  that  are  relevant  to,  and  will  be  used  and  referred  to 
on  a  number  of  occasions  throughout  this  report.  Section  2.1,  presents  a  range  of  algorithms  that 
have  seen  common  use  for  identifying  models  and  features  of  fluids  systems.  Relevant  literature 
and  applications  of  such  methods  are  discussed  in  detail  in  Section  1.2. 

Of  each  of  these  methods,  that  which  is  utilized  more  predominantly  in  this  report  is  DMD, 
which  is  featured  to  various  extents  in  each  of  chapters  3,  4,  5,  and  6.  Additionally,  Galerkin  pro¬ 
jection,  which  utilizes  POD  modes  as  an  efficient  basis  for  approximating  the  governing  equations, 
is  used  in  Chapter  5  as  a  comparison  to  the  EDMD  method  of  nonlinear  system  identification. 
We  include  a  presentation  of  ERA  both  since  it  is  used  in  Section  2.2,  and  to  highlight  the  simi¬ 
larities  with  DMD;  similarities  which  exist  between  numerous  data-driven  linear  modeling/system 
identification  algorithms. 


2.1  Data-driven  modeling  of  fluids  systems 

2.1.1  Proper  orthogonal  decomposition 

The  goal  of  the  proper  orthogonal  decomposition  (POD)  is  to  obtain  a  set  of  empirical  spatial 
modes  that  optimally  represent  a  given  dataset  from  an  energetic  standpoint.  Assume  that  we  can 
decompose  the  dynamics  of  some  system  u(x,t )  (which  could  be  the  time- varying  velocity  field  of 
a  fluid,  say)  by 

OO 

u(x,t )  =  u0(x)  +  y (2-l.l) 
i= 1 

where  u o(x)  is  some  fixed  (often  average)  data,  and  {ui(x)}?!L1  are  a  set  of  orthonormal  basis 
functions  (modes).  POD  takes  these  modes  to  be  those  which  successively  capture  the  most  energy 
of  the  velocity  field.  Each  POD  mode  Ui  satisfies  the  integral 

/  Rij(x,x')ui(x')dx'  =  Xui(x), 

Jn 

where  Rij(x,y)  =  E[ui(x)uj (y)\ ,  with  E  being  the  expectation.  As  indicated  by  Eq.  (5.2.1), 
POD  is  normally  performed  after  first  subtracting  the  mean  (or  perhaps  an  equilibrium  point) 
from  the  data.  This  approach  has  the  advantage  that  uq  satisfies  the  required  non-homogeneous 
boundary  conditions,  meaning  that  all  other  modes  n*  will  satisfy  homogenous  boundary  conditions, 
so  any  linear  combination  of  modes  of  the  form  given  by  Equation  (5.2.1)  will  automatically  satisfy 
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the  correct  boundary  conditions  of  the  problem  at  hand.  In  discrete  terms,  if  we  arrange  finite¬ 
dimensional  data  collected  from  a  simulation  or  experiment  into  a  matrix  Y#,  with  each  column 
representing  a  snapshot  of  data  at  a  given  time,  then  the  POD  modes  are  the  columns  of  U 
in  the  singular  value  decomposition  Y #  =  UYV* .  Here  the  ith  entry  of  the  diagonal  matrix 
Y,  corresponds  to  the  energy  contained  in  the  ith  POD  mode.  In  this  discrete  formulation,  for 
simplicity  we  are  omitting  any  rescaling  of  the  data  that  should  be  performed  so  that  the  modes 
are  orthonormal  with  respect  to  the  usual  inner  product.  That  is,  if  Ui  and  Uj  are  columns  of  U, 
then  we  really  should  have 


u*  ( x')u.i{x')dx ' 


n 

E 

k= 1 


u*(xk)ui(xk)dxk  =  5, 


ij, 


rather  than  u* ut  =  8ij.  The  original  data  Y #  can  then  be  represented  in  terms  of  POD  coefficients 
by  Y #  =  U*Y ft.  If  we  wish  to  reduce  the  dimension  of  this  data,  we  may  do  so  in  an  optimal  way 
(with  respect  to  energy  content)  by  simply  truncating  the  columns  of  U  beyond  a  certain  point, 
which  corresponds  to  removing  POD  modes  that  are  of  sufficiently  low  energy.  Doing  this  gives  a 
reduced  order  approximation  of  the  data  Y#  =  U*Y^ ,  where  Ur  contains  the  first  r  columns  of 
U.  Note  that  there  are  alternative  truncation  techniques  that  may  be  more  effective  than  energy 
maximization  for  certain  applications,  for  example  balanced  POD  [112]  gives  a  reduced  order  linear 
state  space  model  that  is  optimal  with  respect  to  a  given  set  of  sensors  and  actuators. 


2.1.2  Galerkin  projection 

The  idea  behind  GP  is  to  approximate  the  governing  PDEs  that  describe  a  given  system  with  a 
low- dimensional  set  of  ODEs.  This  is  accomplished  by  projecting  the  equations  onto  spatial  POD 
modes  identified  using  the  methods  described  in  Sect.  5.2.1.  We  begin  with  the  incompressible 
Navier-Stokes  equations: 

Du 

— -  =  —  (u  •  V)w  —  uAu  —  VP 

dt  K  ’  (2.1.2) 

V  •  u  =  0. 

If  we  take  the  (spatial)  inner  product  of  Eq.  (5.2.2)  with  a  given  mode  Uj,  we  obtain 

=  —  ((u  ■  V)tt,  Uj)  —  i/  (A u,  Uj)  —  (VP,  Uj) .  (2.1.3) 

Substituting  in  the  finite-dimensional  approximation  of  Eq.  (5.2.1),  we  obtain 

a  =  La  +  Q(a,  a)  +  f,  (2-1.4) 

were  L  is  a  linear  operation  (i.e. ,  a  matrix),  Q  is  a  bilinear  operator  (which  can  be  represented  as 
a  3-tensor),  and  f  is  a  vector,  each  defined  based  on  the  identified  spatial  POD  modes  by 


du 


L%j  —  V  (Auj,  Ui)  ,  Qijk  —  ((^0  '  Hi)  j  f  —  (Vp,  U{)  . 


(2.1.5) 


This  gives  a  means  of  approximating  the  Navier-Stokes  equations  by  a  set  of  nonlinear  ODEs.  As 
mentioned  in  Sect.  6.1,  there  are  many  modifications  that  have  been  proposed  for  this  general 
procedure,  most  typically  to  account  for  the  energy  transfer  to  unmodeled  modes  (i.e.,  the  energy 
cascade  to  finer  spatial  scales).  For  cases  where  spatial  symmetries  exist  (e.g.,  in  the  streamwise 
and  azimuthal  directions  for  circular  pipe  flow),  one  can  show  that  the  POD  modes  must  become 
Fourier  modes,  which  can  simplify  their  computation.  It  is  also  possible  to  “factor  out”  such 
symmetries  by  using  an  optimally  chosen  moving  frame  of  reference  [114,  117]. 
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2.1.3  Dynamic  mode  decomposition 

DMD  has  undergone  a  number  of  formulations,  interpretations  and  modifications  since  its  inception. 
Common  to  all  methods  is  the  requisite  collection  and  arrangement  of  data,  summarized  now. 
Suppose  we  collect  snapshots  of  data  y which  we  assemble  as  columns  in  the  data  matrix  Z. 
For  fluids  systems  y will  typically  be  a  velocity  field  snapshot,  but  more  generally  it  is  a  vector 
of  observations  of  an  evolving  dynamical  system  at  a  given  time.  From  Z ,  we  select  all  pairs  of 
columns  that  are  sampled  at  a  time  difference  At  apart,  and  place  them  into  the  matrices  Y  and 
Y*  (where  the  data  in  a  given  column  of  Y#  was  collected  At  after  the  equivalent  column  of  Y). 
Note  that  if  Z  consists  of  a  sequential  time-series  of  data,  then  Y  and  Y#  are  simply  Z  with 
the  last  and  first  columns  excluded,  respectively.  Let  Y  and  Y#  each  be  n  by  m  matrices,  so  we 
have  m  pairs  of  snapshots,  each  of  size  n.  By  not  explicitly  requiring  a  single  time-series  of  data, 
we  allow  for  larger  or  irregular  time  gaps  between  snapshot  pairs,  the  concatenation  of  data  from 
multiple  time-series,  and  for  the  removal  of  any  corrupted  or  spurious  data.  Recently,  Tu  et  al. 
[149]  proposed  an  interpretation  of  DMD  modes  and  eigenvalues  as  the  eigendecomposition  of  the 
matrix 

A  =  Y#Y+,  (2.1.6) 

where  Y+  denotes  the  Moore-Penrose  pseudoinverse  of  a  matrix  Y.  While  this  is  a  succinct 
interpretation,  and  one  which  will  be  useful  in  the  ensuing  analysis,  it  is  typically  not  an  efficient 
(or  even  feasible)  means  of  performing  DMD  (as  discussed  in  [149]).  This  is  especially  true  when 
n  m.  which  is  often  the  case  in  high-dimensional  fluids  systems.  Instead,  since  Y  and  Y#  have 
rank  at  most  min (m,n),  it  is  typically  more  efficient  to  first  project  the  data  onto  a  subspace  that 
is  (at  most)  of  this  dimension.  One  way  to  do  this  is  by  projecting  the  original  snapshots  onto  the 
POD  modes  of  the  data,  which  is  implicitly  done  in  all  DMD  algorithms.  Note  that  the  POD  modes 
of  Y  are  the  columns  of  U  in  the  singular  value  decomposition  Y  =  UYV*  (though  typically  POD 
is  performed  after  first  subtracting  the  temporal  mean  of  the  data,  which  is  not  done  here).  We 
present  here  a  typical  algorithm  to  compute  DMD,  that  is  most  similar  to  that  proposed  in  [149] 
as  exact  DMD. 

Algorithm  1  (DMD). 

1.  Take  the  reduced  singular  value  decomposition  (SVD)  ofY,  letting  Y  =  UYV* . 

2.  (Optional)  Truncate  the  SVD  by  only  considering  the  first  r  columns  of  U  and  V,  and  the 
first  r  rows  and  columns  of  X  (with  the  singular  values  ordered  by  size),  to  obtain  Ur,  Sr, 
and  Vr 

3.  Let  A  :=  U(Y*VrY-x 

4-  Find  the  eigenvalues  /.q  and  eigenvectors  Wi  of  A,  with  Awi  =  p>iWi, 

5.  Every  nonzero  m  is  a  DMD  eigenvector,  with  a  corresponding  DMD  mode  given  by  ipi  := 
liT1Y#Vr'E-1wi. 

This  method  is  similar  to  the  original  formulation  in  [126],  but  for  the  fact  that  in  step  5 
the  DMD  modes  are  no  longer  restricted  to  lie  within  the  column  space  of  Y.  We  also  explicitly 
provide  the  optional  step  of  truncating  the  SVD  of  Y,  which  might  be  done  if  the  system  is  known 
to  exhibit  low  dimensional  dynamics,  or  in  an  attempt  to  eliminate  POD  modes  that  contain 
only  noise.  We  note  that  this  is  not  the  only  means  to  reduce  the  dimension  of  the  identified 
system  dynamics,  nor  is  it  necessarily  optimal.  Indeed,  [164]  develops  a  method  that  optimizes  the 
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projection  basis  in  parallel  while  performing  a  DMD-like  eigendecomposition.  [71]  takes  a  different 
approach,  seeking  a  small  number  of  nonzero  modes  from  the  full  eigendecomposition  that  best 
approximate  the  system  dynamics.  An  empirical  comparison  between  these  various  dimensionality- 
reduction  techniques  will  be  given  in  Sect  6.3.3.  Note  that  the  continuous  eigenvalues  A  &  of  the 
system  are  related  to  the  discrete  time  eigenvalues  identified  via  DMD  via  A d  =  log(pLi)/ At.  The 
growth  rate  7 j  and  frequency  07  associated  with  DMD  mode  <pi  are  then  given  by  A c*  =  7 ,  +  iui . 

The  matrix  A  is  related  to  A  in  Eq.  (6.2.1)  by  A  =  U*AUr.  While  A  can  be  viewed  as  an 
approximating  linear  propagation  matrix  in  Rn  (i.e.,  the  space  of  original  data  vectors),  A  is  the 
equivalent  propagation  matrix  in  the  space  of  POD  coefficients,  which  we  will  sometimes  refer  to 
as  POD  space.  Another  interpretation  of  A  is  that  it  is  the  spatial  correlation  matrix  between  the 
POD  modes  Ur,  and  the  same  POD  modes  shifted  by  the  assumed  dynamics  A  [126].  If  we  let 
Xfc  =  £/*Xfc  be  the  representation  of  a  given  snapshot  x  in  the  POD  basis  and  let  Y  =  U*Y  and 
Y#  =  U*Y#,  then  it  is  easy  to  verify  that  the  equivalent  of  Eq.  (6.2.1)  in  POD  space  is 

A  =  Y*Y+.  (2.1.7) 

Eq.  (6.2.2)  will  be  useful  for  the  subsequent  analysis  performed  in  this  paper. 

2.1.4  Eigensystem  realization  algorithm 

The  eigensystem  realization  algorithm  (ERA)  is  a  method  that  extracts  a  linear  state  space  model 
from  impulse  response  data.  As  mentioned  in  section  1.2,  it  was  first  conceived  for  analyzing  the 
structural  dynamics  of  spacecraft  in  Juang  and  Pappa  [72],  but  also  shares  close  similarities  with 
a  number  of  previously  proposed  techniques  [e.g. ,  65,  82].  More  details  about  a  range  of  similar 
methods  and  their  potential  applications  can  be  found  in  [151,  152,  109]. 

The  output  of  ERA  is  a  discrete-time  linear  state-space  system,  taking  the  form. 


Xk+i  =  Adxk  +  Bduk 
yk  =  C'rfXk  +  Dduk, 


(2.1.8) 


The  ERA  algorithm  proceeds  as  follows: 

1.  Collect  output  data  from  an  impulse  response  of  the  form  {yo,  yp,  y2P>  •  •  • ,  ymp}  and  {yi,  yp+i,  y2P+i, 

2.  Assemble  the  block  Hankel  matrices 


yo 

yp 

y  2P 

y  mc 

H 

'  = 

yp 

y2P 

y3P 

y  (mc+l) 

1 

_ym0p 

y  (m0+l)P 

y  (m0+2)P 

y  (m0+mc)P_ 

yi 

yp+i 

y2P+i 

ymcp+i 

= 

yp+i 

y2P+i 

y3P+i 

y(mc+l)P+l 

ym, 

oP  +  1 

y  (m0+l)P+l 

y  (m0+2)P +1 

y  (m0+mc)P+l_ 

where  mc  and  m0  are  chosen  such  that  mc  +  m0  <  rn. 
3.  Compute  the  (reduced)  SVD  H  =  UYV1 . 
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4.  Truncate  the  SVD  by  only  considering  the  first  r  columns  of  U  and  V ,  and  the  first  r  rows 
and  columns  of  X  (with  the  singular  values  ordered  by  size),  to  obtain  Ur.  Er,  and  Vr,  where 
r  is  the  desired  model  order. 

5.  The  matrices  of  the  reduced-order  model  of  a  system  with  p  inputs  and  q  outputs  are  given 
by 

Ar  =  S-1/2t^’ffVrSj/2> 

B,  =  the  first  p  columns  of  E}J2Vj , 

I  /2  *  V 

Cr  =  the  first  q  rows  of  UrEr'  , 

A  =  yo- 

Note  that  when  p  =  1,  the  data  pairs  in  step  1  can  just  be  taken  from  an  impulse  response  sequence 
with  its  last  and  first  entries  removed.  This  more  general  formulation  allows  for  the  skipping  of 
data  when  assembling  H ,  which  can  reduce  computational  costs,  while  still  allowing  data  to  be 
used  across  a  large  total  time  window. 

In  general,  input-output  data  might  not  be  available  in  the  form  of  an  impulse  response,  in 
which  case  other  more  general  subspace  methods  may  be  used  (e.g.,  Verhaegen  and  Dewilde  [151]). 
Another  approach  is  to  use  a  technique  such  as  observer /Kalman  filter  identification  [74]  to  com¬ 
pute  an  impulse  response  from  input-output  data  with  random  inputs,  before  applying  ERA.  We 
note  that  there  are  close  similarities  between  DMD  and  ERA.  Indeed,  the  two  algorithms  become 
equivalent  when  mo  =  1,  in  the  sense  that  the  A  matrices  from  either  method  are  only  different  by 
a  similarity  transform  [149]. 


2.2  Motivating  results 


Here,  we  present  some  preliminary  results  that  will  serve  as  motivation  for  the  research  in  the 
following  chapters.  The  main  idea  will  be  to  use  ERA  to  identify  models  for  a  pitching  airfoil 
system,  and  use  these  models  to  design  feedback  controllers  that  allow  for  the  lift  to  be  controlled. 
The  core  methodology  behind  these  results  is  given  in  Brunton  et  al.  [23].  For  brevity,  we  defer  a 
detailed  discussion  of  the  experimental  and  numerical  methods  to  Chapters  3  and  4  respectively. 

We  consider  an  airfoil  undergoing  simple  pitching  motion  about  the  quarter  chord,  with  the 
system  input  being  the  kinematics  of  the  airfoil  (captured  by  the  angular  acceleration,  a)  and  the 
output  the  lift  coefficent,  Cl  =  0  ~  Models  (in  the  form  of  low  order  state-space  realizations 

of  the  system)  are  identified  by  applying  the  eigensystem  realization  algorithm  to  discrete-time 
impulse  response  data.  This  occurs  after  first  extracting  the  components  of  the  lift  that  are  directly 
proportional  to  the  angle  of  attack  and  its  derivatives  (a,  a ,  and  a),  described  in  further  detail  in 
[?  ].  For  direct  numerical  simulations,  such  impulse  response  data  is  directly  simulated,  while  in 
experiments  it  is  acquired  by  applying  the  observer/Kalman  filter  identification  algorithm  (OKID) 
to  the  input /output  data  from  pseudo-random,  frequency  rich  maneuvers.  Feedback  controllers  are 
designed  using  loopshaping,  with  a  desired  loopshape  given  by 


180a2(s  +  1.5a) 
s2(s  +  30a) 


(2.2.1) 


where  the  parameter  a  determines  the  bandwidth  of  the  controller.  This  procedure  can  be  applied 
using  both  theoretical  models  (such  as  the  Theodorsen  model)  and  the  previously  discussed  reduced 
order  models,  identified  in  both  direct  numerical  simulations  (DNS)  and  wind  tunnel  experiments. 
Simulations  are  performed  on  a  two-dimensional  flat  plate  airfoil  at  a  Reynolds  Number  of  100, 
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Figure  2.1:  (left)  Step  responses  of  various  magnitude  for  the  closed-loop  system  in  DNS,  and  (right) 
Normalized  step  responses.  Linear  behavior  in  the  output  (lift  coefficient)  is  observed,  despite  large 
nonlinearities  present  in  the  system,  which  is  evident  from  the  variation  in  the  angle  of  attack  plots. 


using  an  immersed  boundary  projection  method  ([139,  34]).  Experiments  were  conducted  using  a 
NACA0006  airfoil  of  chord  length  0.245m  with  a  free  stream  velocity  of  3  m/s,  giving  a  Reynolds 
number  of  approximately  50,000.  Gusting  conditions  are  generated  in  the  wind  tunnel  by  using  a 
series  of  shutters  downstream  of  the  test  section.  In  spite  of  the  differences  in  parameters  and  condi¬ 
tions  between  the  DNS  and  experimental  work,  we  demonstrate  that  the  same  control  methodology 
can  be  successfully  applied  in  both  cases. 

Feedback  control  is  implemented  for  tracking  reference  lifts  of  a  range  of  magnitudes,  both 
with  and  without  the  presence  of  gusting  disturbances.  To  begin  with,  we  consider  tracking  step 
changes  in  reference  lift  in  DNS.  Figure  2.1  shows  accurate  tracking  of  the  desired  lift  over  a  range 
of  magnitudes,  even  significantly  beyond  the  maximum  value  that  can  be  held  in  steady  conditions 
( Cl  =  0.97).  This  highlights  one  of  the  major  benefits  of  using  feedback  control:  even  though  the 
system  is  clearly  non-linear,  we  are  able  to  make  the  output  (which  is  often  what  we  care  most 
about)  behave  linearly.  Here  the  nonlinear  effects  are  compensated  for  by  the  controller  modifying 
the  input  to  the  system.  Figure  2.2  shows  lift  tracking  step  responses  for  experimental  conditions. 
The  presence  of  noise  and  time  lags  in  the  system  degrades  the  performance.  Nonlinear  effects 
further  limit  the  performance  for  high  amplitude  steps. 

Figure  2.3  shows  experimental  data,  where  feedback  control  is  used  to  track  a  constant,  high 
magnitude  lift  in  the  presence  of  10%  periodic  fluctuations  in  the  freestream  velocity.  Despite 
not  receiving  any  direct  information  about  the  disturbance,  the  controller  is  able  to  correct  for 
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Figure  2.2:  Step  responses  of  various  magnitude  for  the  closed-loop  system  in  wind  tunnel  experi¬ 
ments.  The  controller  performance  degrades  for  high-amplitude  steps 


- measured  controlled 

- predicted  uncontrolled 
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Figure  2.3:  Tracking  a  lift  coefficient  of  0.8  in  wind  tunnel  experiments  in  the  presence  of  periodic 
gusts. 


the  changes  in  lift  coefficient  that  are  caused  by  the  fluctuating  free  stream  velocity,  in  order  to 
maintain  a  nearly  constant  lift.  The  performance  closely  agrees  with  predictions  that  can  be  made 
using  identified  models  for  the  gust  disturbance. 

Having  validating  that  feedback  control  can  be  successfully  implemented  on  simple  step  maneu¬ 
vers,  we  proceed  to  investigate  more  complex  desired  lift  profiles.  Figure  ??  shows  the  performance 
of  the  DNS  system  in  tracking  a  sinusoidally  varying  reference  lift.  Interestingly,  we  note  that  the 
addition  of  a  periodic  component  to  the  reference  lift  at  certain  frequencies  allows  for  successful 
tracking  of  a  higher  average  lift.  In  particular, 

While  these  results  show  one  use  for  reduced  order  models,  they  also  suggest  a  few  limitations 
of  such  an  approach.  Firstly,  accurate  lift  tracking  was  only  possible  when  an  accurate,  real-time 
lift  measurement  was  available.  This  is  primarily  due  to  the  fact  that  the  system  is  nonlinear  . 
This  motivates  the  development  of  nonlinear  modeling  procedures  considered  in  Chapters  3  and  5. 
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Figure  2.4:  Controller  performance  in  tracking  sinusoidally  varying  reference  lift  coefficient. 


Figure  2.4  exposes  some  interesting  flow  physics  for  this  particularly  pitching  airfoil  system, 
where  higher  average  lift  coefficients  are  found  when  pitching  at  a  preferred  frequency.  This  phe¬ 
nomena  is  explored  systematically  in  Chapter  4. 

We  also  observe  a  degradation  of  performance  when  using  experimental,  rather  than  numerical 
data.  Aside  from  issues  with  implementing  feedback  control  in  experimental  systems  that  contain 
time  lags  between  command  and  actuation,  one  cause  of  this  degradation  is  the  presence  of  noise 
in  the  data.  This  is  particularly  important  when  it  is  data  that  is  being  used  to  identify  models. 
We  study  the  presence  of  noise  in  the  dynamic  mode  decomposition  in  Chapter  6,  and  propose 
modified  algorithms  that  can  give  improved  performance  with  noisy  data. 
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Chapter  3 

A  data-driven  modeling  framework 
for  predicting  forces  and  pressures  on 
a  rapidly  pitching  airfoil 


This  work  formulates  a  switched  linear  modeling  procedure  to  understand  and  predict  the  unsteady 
aerodynamic  forces  arising  from  rapid  pitching  motion  of  a  NACA  0012  airfoil  at  a  Reynolds  number 
of  50,000.  The  system  identification  procedure  applies  a  generalized  dynamic  mode  decomposition 
algorithm  to  time-resolved  wind  tunnel  measurements  of  the  lift  and  drag  forces,  as  well  as  the 
pressure  at  six  locations  on  the  suction  surface  of  the  airfoil.  Linear  state  space  models  are  identified 
for  5-degree  pitch-up  and  pitch-down  maneuvers  within  an  overall  angle  of  attack  range  of  0°-20°. 
The  identified  models  accurately  capture  the  effects  of  flow  separation  and  leading-edge  vortex 
formation  and  convection.  It  is  shown  that  switching  between  different  linear  models  can  give 
accurate  prediction  of  the  nonlinear  behavior  that  is  present  in  high- amplitude  maneuvers.  The 
models  are  accurate  for  a  wide  range  of  motions,  including  pitch-and-hold,  sinusoidal,  and  pseudo¬ 
random  pitching  maneuvers.  Providing  the  models  access  to  a  subset  of  the  measured  data  channels 
can  allow  for  improved  estimates  of  the  remaining  states  via  the  use  of  a  Kalman  filter,  which  could 
be  of  use  for  aerodynamic  control  applications. 

3.1  Introduction 

The  flight  of  small,  highly  maneuverable  aircraft,  whether  biological  or  manmade,  is  greatly  im¬ 
pacted  by  unsteady  aerodynamic  effects,  which  can  be  either  beneficial  or  detrimental  to  flight. 
Accurate  understanding  of  such  effects  can  allow  for  the  design  of  aircraft  that  are  more  efficient, 
responsive,  and  robust. 

The  need  to  account  for  unsteady  effects  has  been  recognized  since  soon  after  the  breakthrough  of 
powered  manmade  flight,  in  the  classical  works  of  Wagner  [155],  Theodorsen  [142],  and  Garrick  [50]. 
Indeed,  many  failed  attempts  at  flight  can  probably  be  attributed  to  a  severe  lack  of  understanding 
of  how  to  utilize  such  effects.  These  classical  models  give  significant  insight  into  the  fundamental 
flow  physics  associated  with  unsteady  flight,  such  as  relative  contributions  to  lift  of  the  added 
mass,  quasi-steady  bound  circulation,  and  wake  vortices.  While  such  models  can  be  quantitatively 
accurate  for  cases  of  attached  flow  where  viscous  effects  are  negligible,  they  quickly  lose  validity 
when  dealing  with  separated  flows,  which  are  often  encountered  in  the  extreme  motions  that  are 
possible  for  birds,  insects,  and  micro  and  unmanned  aerial  vehicles  (MAV  and  UAV).  It  is  precisely 
in  these  extreme  cases  that  accurate  predictive  models  are  essential  to  prevent  catastrophic  failure 
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and  ensure  ongoing  successful  flight.  While  more  accurate  predictions  can  be  attained  from  high- 
fidelity  simulations,  the  computational  cost  typically  prohibits  the  direct  use  of  such  simulations 
for  real-time  prediction  and  control. 


Biological  examples  such  as  insects  [18,  120,  156]  and  birds  [153]  have  seemingly  evolved  to  take 
advantage  of  the  high  transient  lift  force  that  can  be  generated  due  to  the  formation  of  a  leading 
edge  vortex  (LEV)  during  rapid  pitch-  up  motion,  for  example.  While  these  give  motivating  ex¬ 
amples  of  the  advantages  of  accurate  understanding  of  unsteady  aerodynamic  effects,  the  preferred 
wing  kinematics  arising  from  evolution  is  highly  specific  and  coupled  to  the  geometry  and  other 
physiological  features  of  the  animal.  Indeed,  the  characteristics  of  unsteady  aerodynamic  effects, 
particularly  for  separated  flows,  seem  to  be  quite  sensitive  to  both  the  geometry[79]  and  Reynolds 
number  [166]  of  the  airfoils.  Studies  into  low  Reynolds  number  flow  over  stationary  [166,  4,  165] 
and  pitching  [154,  1,  131,  30,  80]  symmetric  airfoils  have  revealed,  for  example,  complex  effects 
associated  with  the  stability  and  separation  of  the  suction  surface  boundary  layer,  which  are  again 
highly  sensitive  to  Reynolds  numbers.  These  observations  motivate  the  development  of  general 
modeling  procedures  that  can  be  easily  applied  to  a  range  of  parameter  cases.  In  addition,  it  is 
desirable  for  such  methods  to  be  sufficiently  general  such  that  they  can  be  applied  to  more  realistic 
aircraft  configurations,  rather  than  just  airfoils.  As  an  example,  such  data  driven  modeling  was 
considered  for  the  case  of  accurate  prediction  and  control  of  lift  for  a  low  Reynolds  number  pitching 
airfoil  [22,  23],  using  the  eigensystem  realization  algorithm  [72]  (ERA)  and  observer /Kalman  filter 
identification  [74]  (OKID).  There  has  also  been  a  significant  amount  of  work  in  terms  of  nonlinear 
modeling,  ranging  from  low  order  state-  space  models  formulated  from  theoretical  considerations 
[54],  to  Volterra  series  models  that  have  been  used  to  model  a  range  of  unsteady  aerodynamic  and 
aeroelastic  phenomena[132,  85,  11]. 


More  generally,  rapid  advances  in  both  computational  power  and  experimental  equipment  has 
seen  a  large  increase  in  the  amount  of  data  that  can  be  generated  by  researchers  in  fluid  mechanics. 
This  has  lead  to  the  increased  popularity  of  techniques  such  as  proper  orthogonal  decomposition 
[66]  and  dynamic  mode  decomposition  [126]  (DMD),  which  can  be  useful  to  extract  tractable  models 
and  physical  insight  from  large  fluids  datasets. 


In  the  present  work,  we  use  a  variant  of  DMD  to  identify  linear  state  space  models  for  a  variety  of 
pitch-up  and  pitch-down  maneuvers.  DMD  was  first  introduced  to  the  fluids  community  as  a  means 
to  extract  dynamic  information  from  data  [125],  and  has  subsequently  been  successfully  applied  to 
a  range  of  numerical  and  experimental  fluids  datasets  [126,  128,  127,  70].  Numerous  theoretical 
developments  have  highlighted  DMD’s  connections  to  both  Fourier  [27]  and  Koopman  [119]  spectral 
analyses,  and  other  system  identification  techniques  such  as  ERA  [149].  DMD’s  formulation  has  also 
been  generalized  to  allow  for  the  identification  of  systems  with  inputs  [105],  and  it  is  this  framework 
that  we  make  use  of  in  this  work.  One  advantage  of  the  present  modeling  approach  is  that, 
unlike  those  generated  using  ERA/OKID  the  model  states  can  be  directly  related  to  measurements. 
This  can  allow  easy  switching  between  neighboring  linear  models,  which  subsequently  permits  the 
formation  of  a  switched  linear  model  that  is  capable  of  predicting  nonlinear  behavior.  Our  algorithm 
is  described  in  section  3.2,  which  is  followed  by  a  description  of  the  experimental  setup  in  section 
3.3.  Section  5.4  demonstrates  that  the  obtained  models  are  accurate  for  a  range  of  high-anrplitude 
pitching  maneuvers.  Section  5.5  contains  a  more  general  discussion  of  the  results  and  subsequent 
conclusions  of  this  study. 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


3.2.  SYSTEM  IDENTIFICATION  METHOD 


19 


3.2  System  identification  method 

We  use  a  variant  of  dynamic-mode  decomposition  to  identify  models  describing  the  pitching  airfoil, 
which  is  briefly  described  here.  The  goal  will  be  to  identify  a  family  linear  systems  of  the  form 

Xj-|_i  =  A^yL'i  -f-  -Bfcttj.  (3.2.1) 

By  identifying  different  and  matrices  for  different  angles  of  attack  and  directions  of  pitching, 
we  can  assemble  a  family  of  linear  models  {A^,  Bk}^=l  that  can  subsequently  be  pieced  together  to 
allow  for  accurate  prediction  of  maneuvers  spanning  a  wider  range  of  angles  of  attack  than  any  single 
linear  model  would  be  capable  of.  We  now  describe  in  general  terms  the  identification  procedure 
for  a  linear  model  {A,  B}.  Further  details  concerning  the  specific  data  chosen  to  constitute  the 
state  x  will  be  given  in  Section  3.4.1. 

Suppose  we  collect  a  time-series  of  measurements  Xj,  which  we  assemble  as  columns  into  a 
matrix  X.  From  X .  we  select  all  pairs  of  data  that  are  separated  by  some  nominal  time  At, 
which  we  assemble  into  matrices  X\  and  X2.  If  X  consists  of  uniformly  sampled  data,  then 
X\  and  X2  are  X  with  the  last  and  first  columns  removed,  respectively.  Standard  DMD  can 
be  characterized  as  finding  the  eigendecomposition  of  a  matrix  A  satisfying  (or  approximately 
satisfying)  X2  =  AX i[149].  Depending  on  the  size  of  X\  and  X2,  A  is  either  the  (Frobenius) 
minimum-norm  solution  (if  the  data  matrices  have  more  rows  than  columns),  or  the  least-squares 
solution  (otherwise) .  The  usefulness  and  validity  of  this  approach  relies  assumption  that  the  system 
is  autonomous,  and  not  greatly  affected  by  external  inputs.  If  we  have  known  inputs  Ui  assembled 
into  a  matrix  U,  then  it  is  possible  to  modify  DMD  [105]  to  instead  seek  the  matrices  A  and  B 
satisfying 

X2  =  [A  B]  [Xi  U]  .  (3.2.2) 

Provided  that  the  size  of  the  state  m  is  not  excessive,  we  may  compute  the  augmented  system 
matrices  [A  B]  through 

[AB]=X2[X1U}+ ,  (3.2.3) 

where  +  denotes  the  Moore-Penrose  pseudoinverse.  Since  fluids  systems  are,  in  general,  nonlinear, 
the  ability  of  the  identified  linear  system  in  accurately  modeling  all  data  may  be  limited.  However, 
an  intelligent  selection  of  state  variables  x  can  go  a  considerable  way  towards  factoring  out  much 
of  the  nonlinearity  in  the  system.  To  begin  with,  rather  than  directly  using  force  and  pressure 
measurements,  we  can  instead  consider  deviations  from  the  equilibrium  (or  mean)  values  at  a  given 
angle  of  attack.  This  allows  for  the  resulting  linear  model  to  be  accurate  despite  nonlinear  static 
behavior. 


3.3  Experimental  method 

Experiments  were  conducted  at  the  Andrew  Fejer  wind  tunnel  at  the  Illinois  Institute  of  Technology, 
with  a  diagram  of  the  airfoil  mounting  shown  in  figure  3.1.  A  NACA  0012  airfoil  of  chord  length 
c  =  0.245m  was  used  in  a  test  section  of  length  3m  and  cross-section  0.6m  by  0.6m.  The  airfoil 
spanned  the  width  of  the  test-section,  thus  minimizing  three-dimensional  effects.  The  airfoil  was 
mounted  upon  a  six-axis  ATI  nanol7  force  transducer,  which  allowed  for  the  measurement  of  time- 
resolved  forces  and  moments.  This,  in  turn,  was  mounted  upon  two  pushrods  actuated  by  Copley 
servo  tubes,  allowing  for  pitching  motion  to  be  commanded.  For  the  results  presented  here,  only  the 
rear  pushrod  was  actuated,  which  resulted  in  pitching  about  an  axis  0.11c  from  the  leading  edge. 
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Figure  3.1:  Schematic  of  the  experimental  setup. 


Six  pressure  taps  were  installed  along  the  chord  of  the  airfoil  at  one  spanwise  location,  located  at 
distances  of  0.050c,  0.217c,  0.385c,  0.552c,  0.720c,  and  0.887c  aft  of  the  leading  edge.  The  freestream 
velocity  was  measured  using  a  pitot  tube  and  remained  constant  at  a  nominal  value  of  3m/s,  giving 
a  Reynolds  number  of  approximately  50,000,  and  a  convective  time  tc  =  fj  =  0.0817s.  Note  that 
some  blockage  effects  meant  that  the  freestream  velocity  changed  by  a  small  amount  as  the  angle 
of  attack  changed  (approximately  3%  when  pitching  between  0°  and  20°).  All  forces  and  pressures 
were  nondimensionalized  using  the  averaged  velocity  at  the  relevant  phase  of  the  maneuver.  (Note 
that  this  neglects  unsteady  effects  associated  with  the  changing  velocity,  but  since  the  changes  are 
small,  these  effects  should  be  negligible.) 

Force  and  pressure  data  was  acquired  at  a  frequency  of  1000Hz.  For  each  maneuver,  data  was 
phase-averaged  over  at  least  50  cycles  to  reduce  the  effect  of  measurement  noise.  All  maneuvers 
were  also  performed  with  the  wind  tunnel  off  before  and  after  data  was  collected  with  the  tunnel 
switched  on.  These  results  were  also  phase-averaged,  and  subtracted  from  the  tunnel-on  data. 
This  eliminates  (for  the  force  readings)  the  effects  of  the  mass  of  the  wing,  the  added-mass  terms 
associated  with  accelerating  the  surrounding  air,  and  also  any  other  effects  on  the  measurement 
equipment  resulting  directly  from  the  maneuver  performed.  By  eliminating  added-mass  terms,  we 
isolate  the  circulatory  fluids  forces  arising  from  a  given  pitching  maneuver. 

3.4  Results 

Here  results  are  presented  for  the  identification  (Section  3.4.1)  and  performance  of  the  suite  of 
identified  models.  To  test  the  performance  of  the  family  of  models  that  have  been  identified, 
we  analyze  their  ability  to  predict  a  range  of  other  maneuvers.  These  range  from  compositions 
of  similar  individual  maneuvers  (Section  3.4.2),  to  sinusoidal  (Section  3.4.3)  and  pseudo-random 
(Section  3.4.4)  pitching  maneuvers.  The  latter  two  of  classes  of  maneuver  bear  little  similarity  to 
the  maneuvers  used  for  identification.  In  this  sense,  we  will  be  able  to  show  the  generality  of  these 
models,  which  highlights  that  the  identified  models  represent  more  than  simply  fits  to  the  data, 
and  have  predictive  capabilities. 

3.4.1  System  identification  results 

Models  were  identified  separately  from  pitch- up  and  pitch-down  maneuvers  between  0°-5°,  5°-10°, 
10°-15°,  and  15°-20°,  with  model  states  obtained  from  the  6  pressure  readings  and  the  lift  and 
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drag  measurements.  The  prescribed  maneuvers  take  the  canonical  form[47] 

MG  cosh(a(t  —  ti))  cosh(at2) 

rnax(G)  ’  ^  cosh(a(t  —  ^2))  cosh(ati) 


(3.4.1) 


Nondimensionalizing  time  by  we  take  t2  —  t\  =  2,  and  a  =  =  122.4  (and  M  =  5°).  With 

these  parameter  values,  a  becomes  the  main  governing  parameter  that  determines  the  rate  of  the 
step. 

The  duration  of  the  pitch  was  approximately  4  convective  times  (4^).  Static  data  at  the  corre¬ 
sponding  angles  of  attack  was  first  subtracted  from  all  measurements,  and  all  data  was  nondimen- 
sionalized  (forces  by  \pcU2  and  pressures  by  \pU2).  To  identify  models,  the  DMD-type  algorithm 
described  above  was  used,  which  allows  for  data  with  inputs  (which  in  this  case  was  taken  to  be 
measurements  of  either  a  and  a,  or  just  a).  Using  this  method,  we  arrive  at  models  of  the  form 


xi+i  =  Ax*  +  But, 

as  described  in  Section  3.2.  Here  we  let  x  =  [Cpi  Cp 2  •  •  •  Cp q  Ci  Cd]T ,  and  u  =  [a  a]T .  Note  that 
while  we  treat  a  as  an  input  for  convenience,  the  fact  that  it  is  entirely  dependent  on  a  means 
that  we  could  also  treat  it  as  an  additional  system  state.  Here  ~  represents  the  deviation  from  an 
equilibrium  condition,  Ci  =  Ci  —  Cf(a),  where  Cf(a )  is  the  equilibrium  value  at  a  given  angle  of 
attack.  We  first  attain  this  equilibrium  data  for  angles  of  attack  in  the  range  a  G  [0°,22°],  which 
is  shown  in  Figure  3.2.  To  motivate  the  development  of  unsteady  models,  we  also  show  how  data 
acquired  for  a  pitching  airfoil  deviates  from  these  equilibrium  values.  Considering  just  the  static 
data,  we  observe  that  the  lift  coefficient  increases  close  to  linearly  (with  a  slope  of  approximately 
1.77 r)  between  0°  and  8°,  before  the  lift  curve  reaches  its  peak  and  then  plateaus  between  10°  and 
15°,  before  again  increasing  beyond  15°.  The  lift  plateau  corresponds  to  the  airfoil  stalling,  with  the 
flow  over  the  suction  surface  becoming  separated.  Further  evidence  for  this  comes  from  examining 
both  the  drag  curve,  which  sees  a  large  increase  in  drag  beyond  a  =  8°,  and  in  the  first  two  pressure 
coefficients,  which  give  a  sharp  drop  in  pressure  beyond  this  angle.  Prior  to  full  separation,  there 
is  evidence  for  partial  separation  towards  the  rear  of  the  airfoil.  Pressure  sensors  3-6  all  measure 
a  drop  in  pressure  at  a  critical  angle  between  2°  and  6°,  which  appears  to  signify  the  separation 
point  moving  upstream  of  the  given  sensor. 

Returning  now  to  the  system  identification  procedure,  Figures  3.3  and  3.4  show  the  performance 
of  each  model  in  predicting  the  pressure  and  force  coefficients  for  the  maneuver  upon  which  they 
were  identified.  For  reference,  the  static  pressure  and  force  coefficients  at  the  instantaneous  angle 
of  attack  are  also  shown  in  Figures  3.3  and  3.4.  Rather  than  subtracting  the  full  static  curves 
before  system  identification,  we  found  improved  results  by  assuming  linear  variation  in  the  static 
values  throughout  the  maneuver.  This  avoids  issues  with  separation-related  “jumps”  occurring  at 
different  angles  of  attack  in  for  the  static  and  moving  airfoil,  which  makes  the  static-subtracted 
data  less  smooth. 

We  can  identify  two  dominant  features  of  the  pitch-up  and  pitch-down  behavior.  A  temporary 
rise  in  Ci,Cd  and  —  Cp  is  observed  which  is  consistent  with  the  formation  and  convection  of  a 
leading-edge  vortex  (see  e.g.  all  measurements  for  pitching  between  15°  and  20°),  and  a  time-lag 
in  reaching  the  steady  state  value,  most  likely  due  to  the  boundary  layer  requiring  time  to  reach 
its  new  equilibrium  configuration  (see  e.g.  Cpq  when  pitching  between  0°  and  5°). 

We  finally  note  that  we  obtain  quite  different  models  for  pitch-up  and  pitch-down  maneuvers. 
To  show  this  explicitly,  Figure  3.5  shows  the  inaccuracy  of  the  prediction  of  a  5-10°  pitch-up 
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Figure  3.2:  Force  and  pressure  coefficient  data  for  a  static  airfoil  at  angles  of  attack  between  0°  and 
22°,  as  well  as  for  an  airfoil  sinusoidally  pitching  between  5°  and  15°  at  a  rate  k  =  =  0.051. 
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Figure  3.3:  Force  and  pressure  data  for  5°  pitch-up  maneuvers  for  starting  angles  of  0°,  5°,  10° 
and  15°,  which  were  used  for  system  identification.  In  all  cases,  the  identified  models  accurately 
replicate  the  experimental  data.  Also  shown  is  the  static  data  at  the  relevant  instantaneous  angle 
of  attack. 
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Figure  3.4:  Force  and  pressure  data  for  5°  pitch-down  maneuvers  for  starting  angles  of  5°,  10°,  15° 
and  20°,  which  were  used  for  system  identification.  In  all  cases,  the  identified  models  accurately 
replicate  the  experimental  data.  Also  shown  is  the  static  data  at  the  relevant  instantaneous  angle 
of  attack. 
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Figure  3.5:  Example  of  inaccurate  prediction  of  a  pitch-up  model  (5  —  10°  model,  red  curves)  on 
pitch  down  maneuver  (10  —  5°  maneuver,  blue  curves). 


model  for  a  10-5°  pitch-down  maneuver,  which  arises  primarily  because  the  5-10°  pitch-up  model 
predicts  the  existence  of  a  time-delay,  which  is  not  present  in  the  10-5°  pitch-down  data.  This 
has  important  implications  for  the  use  of  pseudo-random  system  identification  maneuvers,  which 
necessarily  incorporate  both  pitching  up  and  pitching  down  motion. 

3.4.2  Multiple  pitch-up  and  pitch-down  maneuvers 

We  now  consider  a  maneuver  consisting  of  two  pitch-ups  followed  by  two  pitch-downs,  each  in  rapid 
succession.  We  attempt  to  predict  the  maneuver  by  switching  between  the  relevant  models  for  each 
pitch-up  and  pitch-down.  For  this  maneuver,  four  different  models  are  used.  Given  that  the  state 
of  each  model  consists  of  the  same  variables,  this  is  simply  a  matter  of  switching  the  A  and  B 
matrices  used  to  propagate  the  system. 

The  results  for  this  procedure  in  predicting  Ci,  Cd ,  as  well  as  two  of  the  pressure  coefficients 
are  shown  in  Figure  3.6,  where  we  have  considered  double  pitch-up/down  maneuvers  between  10° 
and  20°  with  different  pitching  rates.  We  vary  the  a  parameter  from  equation  3.4.1:  to  modify  the 
pitching  rate  (halving  and  doubling  it  from  the  value  used  in  system  identification).  In  all  cases, 
we  switch  between  sub-models  at  t  =  10,  20  and  30  convective  times,  using  the  final  predicted  state 
from  one  sub-model  as  the  initial  condition  for  the  next.  To  give  some  basis  for  comparison,  we 
show  the  performance  of  a  single  linear  model  (that  identified  from  a  5-10°  pitch-up)  in  predicting 
this  maneuver  in  Figure  3.7.  We  note  that  the  only  section  of  this  maneuver  that  this  model 
accurately  predicts  is  that  which  is  most  similar  to  its  identification  maneuver. 

Figure  3.8  shows  a  quadruple  pitch-up  and  -down  maneuver,  which  switches  between  all  models. 
From  all  of  these  results,  we  find  that  switching  between  models  generally  works  well,  though 
sometimes  it  can  induce  “jumps”  immediately  after  switching,  particularly  when  switching  between 
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Figure  3.6:  Switched  model  performance  in  predicting  pressure  and  force  coefficients  for  double 
pitch  up/down  maneuvers  with  different  pitch-rates,  between  10°  and  20°.  Switching  between 
sub-models  occurs  at  t  =  10,  20  and  30  convective  times,  using  the  final  predicted  state  from  one 
sub-model  as  the  initial  condition  for  the  next.  The  middle  plot  uses  the  same  pitch-  rate  as  the 
maneuvers  used  for  system  identification. 
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Figure  3.7:  Performance  of  a  single  (5-10°  pitch-up)  model  in  predicting  pressure  and  force  coeffi¬ 
cients  for  double  pitch  up/down  maneuver. 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


3.4.  RESULTS 


27 


100 


100 


Q 

O 


Figure  3.8:  Performance  of  switched  model  in  predicting  pressure  and  force  coefficients  for  quadruple 
pitch  up/down  maneuver.  Switching  between  models  occurs  at  multiples  of  10  convective  times 
between  t  =  10  and  t  =  70,  where  the  final  predicted  state  from  one  sub-model  as  the  initial 
condition  for  the  next. 


the  pitch-up  and  pitch-down  models  between  5°  and  10°.  It  is  possible  that  these  could  be  eliminated 
or  reduced  by  further  refining  the  system  identification  and/  or  switching  procedure. 

3.4.3  Sinusoidal  pitching 

Next,  we  consider  high-amplitude  sinusoidal  pitching  maneuvers,  pitching  between  0°  and  20°  at 
rates  /  =  0.2  Hz  and  0.4  Hz,  giving  a  reduced  frequencies  k  =  ^  =  0.051  and  0.103.  In  Figure 

3.9  we  show  the  predicted  pressures  and  forces  when  using  a  single  model  (arbitrarily  taken  to  be 
the  pitch-up  model  from  5°  to  10°),  a  switched  model,  and  a  switched  equipped  with  a  Kalman 
filter  (that  gives  access  to  the  1st  and  6th  pressure  measurements).  Details  concerning  the  design  of 
this  Kalman  filter  are  given  in  Appendix  1.  We  find  that  the  switched  model  performs  better  than 
any  single  linear  model,  and  that  improved  accuracy  in  all  measurements  can  be  achieved  when 
using  the  Kalman  filter.  The  latter  observation  demonstrates  that,  even  if  the  models  themselves 
have  some  inaccuracies  in  predicting  the  outputs,  access  to  measurements  of  a  subset  of  these 
outputs  can  improve  the  prediction  of  all  outputs.  This  is  relevant  for  the  use  of  such  models  for 
real-time  control,  where,  for  example,  we  may  seek  to  attain  a  desired  lift  force  using  only  pressure 
measurements. 

3.4.4  Pseudo-random  pitching 

We  finally  consider  the  case  where  the  angle  of  attack  varies  in  a  pseudo-random  manner.  Figure 

3.10  shows  the  performance  of  a  switched  model  equipped  with  a  Kalman  filter  in  predicting  the 
pressures  and  forces  for  a  pseudo-random  pitching  maneuver.  Again,  we  observe  close  agreement 
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- Experimental  data 

- Single  model,  5-10  deg 

- Switched  model 

- Switched  model  with  Kalman  filter 
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Figure  3.9:  Predicted  and  actual  pressure  and  force  coefficients  for  high-amplitude  sinusoidal  pitch¬ 
ing  at  dimensionless  frequencies  k  =  ^7^  =  0.051  and  0.103.  Predictions  are  made  using  a  model 
both  with  and  without  a  Kalman  filter.  When  a  Kalman  filter  is  used,  the  model  is  given  access  to 
the  1st  and  6th  pressure  measurements. 


between  the  measured  and  predicted  results.  This  close  agreement  highlights  the  full  generality  of 
the  switched  model,  as  it  is  capable  of  accurately  predicting  the  behavior  of  the  airfoil  forces  and 
pressures  for  arbitrary  high-amplitude  pitching  motions. 


3.5  Discussion  and  conclusions 

The  results  presented  in  Section  5.4  demonstrate  that  the  system  identification  technique  described 
in  Section  3.2  can  be  of  use  for  unsteady  aerodynamic  modeling  applications.  The  fact  that  accurate 
models  were  attained  from  very  simple  pitch-up  and  pitch-down  maneuvers  gives  the  procedure 
an  advantage  over  the  OKID  algorithm,  which  typically  requires  a  concatenation  of  a  variety  of 
motions  to  obtain  accurate  models [22],  The  absence  of  internal  states  in  the  resulting  models  mean 
that  they  are  naturally  suited  for  piecing  together  for  the  formation  of  a  global  switched  model. 
This  process  is  difficult  for  ERA  models,  where  the  internal  states  are  not  directly  associated  with 
physical  measurements.  Having  measurements  directly  associated  with  model  states  means  that  the 
dimension  of  the  observables  must  be  at  least  as  large  as  the  dimension  of  the  underlying  dynamics 
(or  their  approximating  model),  though  this  restriction  could  be  relaxed  if  we  were  to  concatenate 
the  data  with  time-shifted  measurements  (as  is  done  in  ERA),  or  by  using  transformations  of  the 
original  data  [160].  Conversely,  the  fact  that  the  models  are  accurate  suggests  that  8th  order  linear 
models  are  sufficient  to  capture  the  phenomena  present  in  the  maneuvers  considered.  Indeed,  in 
many  cases  it  was  found  that  it  was  possible  to  apply  balanced  truncation  to  reduce  the  dimension 
of  the  identified  models  without  significant  degradation  of  predictive  accuracy. 

In  general,  linear  modeling  techniques  are  appealing  due  to  the  simplicity  of  their  identification 
and  formulation,  and  the  ease  of  use  in  simulation  and  controller  design.  Their  accuracy  in  the 
prediction  of  nonlinear  dynamics,  however,  will  typically  be  fundamentally  limited  to  a  region  in 
phase-space  that  is  near  to  the  identification  maneuver.  Scheduling  between  a  family  of  linear 
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Figure  3.10:  Performance  of  the  switched  model,  equipped  with  a  Kalman  filter  and  measurements 
of  the  first  and  last  pressure  coefficient,  in  predicting  pressure  and  force  coefficients  for  high- 
amplitude  pseudo-random  pitching  maneuver. 
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models  can  go  some  way  to  incorporating  nonlinear  effects  into  a  global  model,  thus  increasing 
the  region  in  phase-space  where  such  models  are  accurate.  This  work  demonstrated  that,  in  the 
a  direction  of  phase-space,  such  an  approach  can  work  between  0°  and  20°,  which  includes  the 
regimes  where  the  flow  over  the  suction  surface  is  completely  attached,  partially  separated,  and 
fully  separated.  This  range  was  the  maximum  available  given  physical  limitations  of  the  airfoil 
mounting,  but  we  imagine  that  separated  flows  at  higher  angles  of  attack  should  also  be  able  to 
be  accurately  modeled,  given  that  they  are  phenomenologically  similar  (in  terms  of  being  fully 
separated)  to  those  near  20°. 

In  terms  of  the  applicability  of  the  model  for  maneuvers  with  different  pitch-rates  (i.e.  the 
a  direction  of  phase  space),  we  see  in  Figure  3.6  that  the  models  remain  accurate  for  a  range  of 
values  of  pitch  rates.  Looking  at  Figure  3.9,  however,  we  see  that  switching  between  models  while 
pitching  at  a  relatively  fast  rate  can  lead  to  some  degradation  of  model  accuracy.  This  is  a  known 
limitation  of  gain  scheduling  models  in  general[130]. 

The  fact  that  the  acquired  data  was  phase-averaged  over  a  number  of  cycles  means  that  any 
unsteady  phenomena  that  are  not  phase-locked  with  the  pitching  motions  will  be  averaged  out  of 
the  identified  models.  Particularly  for  separated  boundary  layers,  such  effects  (which  can  occur 
on  a  faster  timescale  to  the  pitching  motions)  can  be  significant,  even  if  they  are  not  directly 
controllable  by  pitching  motion.  Further  work  could,  for  example,  incorporate  such  dynamics  into 
state  estimators,  which  could  improve  the  real-time  predictive  power  of  such  models. 

The  data  that  is  obtained  for  the  cases  of  a  pitching  and  stationary  airfoil  is  also  of  fundamen¬ 
tal  fluid  mechanical  interest,  which  will  be  further  investigated  by  investigating  the  time- varying 
velocity  field  using  particle  image  velocimetry  (PIV).  Specifically,  it  would  be  interesting  to  explore 
whether  a  small  number  of  measurements  could  be  used  to  accurately  predict  not  only  the  pressures 
and  forces  (as  was  done  in  the  present  work),  but  also  the  entire  velocity  Held  in  the  vicinity  of  the 
airfoil. 


Appendix  1:  Kalman  filter  design 

The  linear  models  (taking  the  form  given  in  equation  3.2.1)  identified  in  this  work  give  a  prediction 
of  forces  and  pressures  given  knowledge  of  the  airfoil  pitching  kinematics.  In  some  situations, 
it  might  be  possible  to  supplement  knowledge  of  these  kinematics  (i.e.,  the  system  inputs)  with 
some  number  of  additional  measurements.  This  section  briefly  describes  the  setup  and  design  of 
a  Kalman  filter  [78],  that  is  used  to  improve  the  estimate  of  the  state  of  the  system  using  such 
additional  measurements.  Suppose  we  have  a  state  space  system  of  the  form 


Xj+i  =  Ax,  +  Bui  +  Gw, 

Yi  =  Cxj  +  Dui  +  HwL  +  Vi, 

which  is  a  generalization  of  equation  3.2.1  to  include  the  influence  of  plant  disturbances  w,  as  well 
as  an  output  equation  that  includes  sensor  noise  v.  For  the  purposes  of  this  work,  we  will  assume 
that  H  =  0,  D  =  0,  and  C  is  a  matrix  that  selects  a  subset  of  the  states  x  as  outputs.  The 
Kalman  filter  gives  an  estimate  of  the  state  x  from  the  system  inputs  u  and  (potentially  noisy) 
measurements  y,  whose  dynamics  are  governed  by  the  equation 

Xj+i  =  Ax,  +  Bui  +  L(yi  -  Cxi  -  Dut ). 

Here  the  matrix  L  gives  the  optimal  state  estimate  for  given  disturbance  and  noise  covariance 
matrices,  Q  =  E( wwT)  and  R  =  £'(vvi),  respectively.  Further  details  concerning  Kalman  filter 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


3.5.  DISCUSSION  AND  CONCLUSIONS 


31 


design  and  the  computation  of  L  may  be  found  in  standard  optimal  control  textbooks[136].  For 
this  work,  we  take  Q  and  R  to  be  appropriately  sized  diagonal  matrices,  and  set  all  diagonal  entries 
to  be  equal  aside  from  the  entries  of  Q  corresponding  to  the  lift  and  drag  states,  which  we  decrease 
by  a  factor  of  10  to  avoid  excessive  oscillations  in  the  estimated  force  coefficients.  We  find  that 
Kalman  filter  performance  is  relatively  insensitive  to  changes  in  these  weights. 
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Chapter  4 


Lift  enhancement  of  high  angle  of 
attack  airfoils  using  periodic  pitching 


In  this  work,  we  study  a  sinusoidally  pitching,  two-dimensional  flat  plate  airfoil  at  a  Reynolds 
number  of  100,  across  a  range  of  pitching  amplitudes,  frequencies,  mean  angles  of  attack,  and  pitch 
axis  locations.  We  report  on  the  lift,  drag,  and  wake  structures  present  in  different  regions  of 
parameter  space.  We  examine  the  average  and  spectral  properties  of  the  forces  on  the  airfoil,  and 
use  dynamic  mode  decomposition  to  examine  the  structures  and  frequency  content  of  the  wake. 
We  give  focus  to  a  number  of  regions  in  parameter  space  where  interesting  behavior  is  observed.  In 
particular,  we  find  that  in  the  regime  where  the  flow  on  the  upper  surface  of  the  airfoil  is  separated, 
but  the  steady  wake  is  stable,  pitching  at  a  specific  frequency  excites  a  vortex  shedding  mode  in 
the  wake,  leading  to  substantial  increase  in  the  lift  and  drag  forces.  This  phenomena  is  insensitive 
to  pitch-axis  location  and  amplitude.  At  higher  angles  of  attack  where  the  wake  for  a  steady  airfoil 
exhibits  periodic  vortex  shedding,  we  find  that,  in  addition  to  this  mean  lift  maxima,  the  interaction 
between  the  natural  and  forced  modes  gives  rise  to  more  complex  behavior. 


4.1  Introduction 

The  unsteady  motion  of  airfoils  at  low  Reynolds  number  and  high  angle  of  attack  leads  to  a  range  of 
phenomena  that  cannot  adequately  be  explained  by  classical  aerodynamic  theories.  It  is  precisely 
these  conditions  that  are  encountered  by  small  fliers,  ranging  from  biological  examples  such  as 
birds [153],  insects[18,  120,  156],  and  bats,  or  manmade  UAVs  and  MAVs.  The  vortex  dynamics 
excited  by  airfoil  motion,  actuation,  or  indeed  present  in  the  natural  flow  at  sufficiently  high  angles 
of  attack,  can  significantly  affect  aerodynamic  performance  [87].  This  motivates  work  that  seeks 
to  understand  and  control  such  phenomena,  and  indeed  modeling  the  dynamics  of  pitching  and 
plunging  airfoils  has  attracted  significant  recent  attention[22,  23,  10,  57,  61,  39,  168]. 

This  work  will  investigate  the  interaction  between  periodic  vortex  shedding  that  can  occur  for 
bluff  bodies  (or  airfoils  at  sufficiently  large  angle  of  attack),  and  imposed  sinusoidal  pitching  motion. 
In  the  case  of  plunging  motion,  it  has  been  observed  that  lock  on[121]  can  occur  between  plunging 
frequency  and  natural  vortex  shedding[167,  31,  33],  while  similar  phenomena  have  been  found  for 
surging  oscillations  over  a  wide  range  of  Reynolds  number  [29].  We  will  show  that  similar  phenomena 
are  observed  in  the  case  of  pitching  motion.  In  experimental  conditions,  plunging  oscillations  may 
also  lead  to  a  bifurcation  of  the  wake  direction  [32], 

It  has  been  observed  that  optimal  pitching  trajectories  maximize  the  circulation  that  is  entrained 
in  leading  edge  vortices[90],  linking  the  concepts  explored  here  to  the  general  notion  of  a  formation 
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Figure  4.1:  Computational  domain  used  for  this  study,  with  the  size  and  location  of  the  airfoil 
shown.  Dashed  lines  represent  the  borders  of  each  nested  grid. 


number  that  can  be  used  to  explain  characteristic  sizes  and  frequencies  associated  with  vortex 
phenomena[51]. 

While  this  work  will  only  consider  two-dimensional  airfoils,  we  note  that  three  dimensional 
phenomena  can  give  rise  to  additional  complexity.  For  example,  aspect  ratio  effects  being  significant 
on  lift  enhancement  due  to  periodic  forcing[140],  and  in  the  case  of  three  dimensional  airfoils  that 
are  free  to  exhibit  rolling  motion,  forced  pitching  can  lead  to  self-excited  roll  oscillations [144]. 

There  are  strong  parallels  between  studying  the  effect  on  lift  and  drag  in  the  context  of  lifting 
airfoils,  and  in  the  investigation  of  thrust-generating  airfoils [42] ,  where  it  is  found  that  propulsive 
efficiency  is  maximized  when  flapping  produces  a  reverse  von-Karman  wake  that  excites  the  least 
stable  spatial  mode  of  the  mean  wake  flow[146,  145].  We  describe  the  numerical  method  and  scope 
of  the  work  in  section  4.2,  before  presenting  our  results  in  section  5.4.  We  will  focus  on  presenting 
and  analyzing  results  at  parameters  where  the  pitching  motion  triggers,  strengthens,  and  interferes 
with  vortex  shedding. 

4.2  Numerical  method  and  scope  of  investigation 

We  use  an  immersed  boundary  projection  method[139,  34]  to  perform  direct  numerical  simulations 
of  the  incompressible  Navier-Stokes  equations.  The  domain  consists  of  four  nested  grids  about  a 
flat  plate  airfoil.  A  diagram  of  the  computational  domain  is  shown  in  Figure  5.1.  Each  of  the 
four  grids  contains  600  by  300  grid  points,  with  a  total  computational  domain  extending  96  and 
48  chord  lengths  in  the  streamwise  and  transverse  directions,  respectively.  The  Reynolds  number 
(based  on  chord  length  and  freestream  velocity)  is  fixed  at  100  throughout.  The  relative  compu¬ 
tational  cheapness  of  such  two-dimensional,  low  Reynolds  number  configurations  makes  thorough 
investigations  of  high-dimensional  parameter  spaces  feasible.  Resolution  studies  were  performed  to 
ensure  that  the  resolution  and  extent  of  the  domain  were  sufficient.  To  give  an  example  of  this, 
we  show  in  Figure  ??  results  from  applying  Crank-Nicholson  and  third-order  Runge-Kutta  time 
steppers  were  used  to  evolve  the  linear  and  nonlinear  terms  respectively,  with  step  size  ranging 
between  At  =  0.0005 c/U  and  At  =  0.01  c/U,  with  the  smaller  range  of  At  required  to  resolve 
pitching  motions  with  larger  frequencies  and/or  amplitudes.  In  all  cases,  we  run  the  simulations 
for  sufficiently  long  such  that  any  limit  cycle  or  long-time  behavior  is  reached  before  the  data  that 
is  used  for  analysis  is  collected.  We  consider  airfoil  kinematics  of  the  form 

a(t)  =  ckm  +  ola  sin(2-7r/*t),  (4.2.1) 
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Figure  4.2:  Resolution  study  to  determine  appropriate  size  of  computational  domain 


where  f*  =  fc/U  is  a  dimensionless  frequency.  We  perform  simulations  with  the  mean  angle  of 
attack  varying  in  5°  increments  between  15°  and  45°,  with  pitching  amplitudes  of  1°,  2°,  5°  and 
10°.  We  consider  frequencies  in  the  range  /*  6  [0.01,5],  with  approximately  20  frequencies  used 
for  each  olm  and  ola-  For  some  cases,  additional  frequencies  are  added  to  improve  local  resolution 
in  parameter  space.  Performing  these  simulations  with  for  pitching  about  the  leading  edge  and 
midchord,  this  results  in  a  total  of  approximately  1120  individual  simulations. 


4.3  Results 

4.3.1  Static  data 

To  give  a  sense  for  the  behavior  of  the  stationary  airfoil,  we  show  a  typical  lift  curve  in  Figure  4.3. 
We  observe  the  expected  linear  relationship  between  the  angle  of  attack,  a ,  and  the  lift  coefficient, 
Cl  for  low  angles  of  attack.  Once  the  angle  of  attack  becomes  sufficiently  large  ( a  >  10°),  flow 
separation  on  the  upper  surface  of  the  airfoil  leads  to  a  shallow  lift  slope.  At  a  ~  20°,  the  flow 
is  separated  and  steady.  Beyond  a  critical  angle  of  attack  ac  ~  27°,  the  steady  solution  becomes 
unstable,  and  periodic  vortex  shedding  is  observed.  This  is  an  example  of  a  supercritical  Hopf 
bifurcation  that  is  seen  in  the  wake  of  bluff  bodies  as  the  Reynolds  number  is  increased[135].  Note 
that  in  the  case  of  an  airfoil  at  an  angle  of  attack,  the  projected  area  csin(a)  is  the  effective  length 
parameter  for  determining  the  location  of  the  bifurcation.  For  a  >  ac,  the  system  exhibits  higher 
lift  (and  drag)  than  would  occur  at  the  unstable  equilibrium  solution. 

4.3.2  Force  analysis 

In  this  section,  we  study  the  lift  and  drag  forces  for  pitching  motion  with  various  amplitudes, 
frequencies,  and  mean  angles  of  attack. 
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Figure  4.3:  Lift  curve  for  stationary  airfoil,  showing  regions  where  the  equilibrium  is  stable  ( a  < 
27°)  and  unstable  ( a  >  27°),  above  which  periodic  vortex  shedding  occurs. 


Figure  4.4  shows  the  mean  lift  coefficient  as  a  function  of  the  dimensionless  frequency,  /*,  for 
pitching  with  a  range  of  amplitudes,  q:a,  and  mean  angles  of  attack,  o>m ■  We  observe  a  distinct 
local  peak  in  Cl  for  all  values  of  am  and  a  a,  aside  from  pitching  with  low  amplitudes  (a  a  =  1°  or 
2°)  about  an  =  15°.  The  location  of  this  lift  peak  moves  slightly  as  am  varies,  from  approximately 
/*  =  0.3  at  Qjf  =  20°,  to  /*  =  0.23  at  =  45°.  For  ocm  >  30°,  we  observe  a  second  peak  emerging 
at  approximately  twice  the  frequency  of  the  dominant  peak.  This  suggests  that  we  excite  dynamics 
that  give  enhanced  lift  when  pitching  at  both  a  fundamental  frequency  and  its  first  harmonic.  We 
note  also  that  the  size  of  the  lift  increment  seems  to  be  largest  for  the  intermediate  base  angles  of 
25°  and  30°. 

Figure  4.5  shows  the  mean  drag  coefficient  for  the  same  range  of  parameters  as  Figure  4.4.  We 
find  that  there  are  increases  in  drag  that  show  similar  behavior  to  those  for  the  lift.  To  compare  the 
changes  in  lift  and  drag  more  explicitly,  we  plot  the  ratio  between  mean  lift  and  drag  coefficients  in 
Figure  4.6.  For  larger  pitching  amplitudes,  there  is  an  increase  in  the  lift-to-drag  ratio  at  certain 
frequencies.  The  frequency  of  the  lift-to-drag  peak  shows  similar  behavior  to  the  peaks  in  both  lift 
and  drag  (note  that,  unlike  the  lift  and  drag  forces,  lift  over  drag  decreases  with  increasing  ajf). 
This  finding  is  potentially  important  for  the  effectiveness  of  such  motions  in  flight. 

Turning  our  attention  back  to  the  lift  coefficient,  we  plot  in  Figure  4.7  the  increase  in  lift 
coefficient  from  the  fixed-wing  value  at  am,  normalized  by  the  amplitude  of  pitching.  For  angles  of 
attack  below  the  critical  angle  (ac)  at  which  vortex  shedding  occurs,  the  lift  increment  is  slightly 
larger  for  larger  a  a,  even  after  normalizing,  despite  the  fact  that  the  mean  lift  is  lower  than  the 
fixed-wing  value  across  other  frequencies.  This  could  be  due  to  the  fact  that,  for  larger  a. a-,  the 
maximum  angle  of  attack  attained  in  a  pitching  cycle  is  closer  to  or  exceeds  ac,  and  thus  better 
able  to  excite  vortex  shedding.  Conversely,  for  am  >  25°,  the  normalized  lift  increment  is  slightly 
larger  for  smaller  pitching  amplitudes,  though  in  some  cases  this  is  in  line  with  the  larger  average 
lift  increments  present  for  lower  pitching  amplitudes  across  all  frequencies. 

To  analyze  the  effect  of  pitch  axis  location,  we  show  in  Figure  4.8  the  normalized  lift  increment 
(i.e.,  the  same  quantity  plotted  in  Figure  4.7)  for  pitching  about  the  leading  edge,  rather  than  the 
midchord.  We  find  very  similar  results  to  pitching  about  the  midchord,  with  maximum  increased 
lift  at  dimensionless  frequencies  between  0.25  and  0.3.  The  lift  increment  is  larger  for  pitching 
about  the  leading  edge,  which  could  be  due  to  the  increased  range  of  motion  of  the  trailing  edge 
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Figure  4.4:  Mean  lift  coefficient  for  a  range  of  pitching  amplitudes  cca,  frequencies,  /*,  and  mean 
angles  of  attack,  a«.  Pitching  is  about  the  midchord. 


for  the  same  angular  pitch  amplitude.  Note  also  that  leading  edge  pitching  will  result  in  larger 
added  mass  forces,  which  for  nonzero  cum  will  increase  lift  and  decrease  drag,  particularly  for  high 
pitching  frequencies.  The  remainder  of  this  section  will  consider  pitching  about  the  midchord.  The 
effect  of  pitch  axis  location  will  be  investigated  more  thoroughly  in  the  next  section. 

Figure  4.9  indicates  how  the  mean  lift  compares  to  the  maximum  and  minimum  lift  for  a  a  =  5°. 
Note  that  in  some  of  the  cases  the  lift  is  not  periodic  with  the  period  of  forcing,  so  these  maximum 
and  minimum  values  are  the  global  extrema  over  many  cycles.  We  observe  (particularly  clearly 
for  lower  mean  angles  of  attack)  that  as  the  frequency  increases  from  low  values,  the  amplitude  of 
the  lift  response  increases,  while  the  mean  remains  approximately  constant.  Above  f*  ~  0.2,  the 
amplitude  of  the  variation  in  lift  decreases,  but  with  the  lift  minimum  rising  more  abruptly  than 
the  lift  maximum  falls.  This  asymmetry  produces  the  higher  average  lift  that  is  observed  in  the 
range  0.1  <  f*  <  0.5.  For  higher  frequencies,  larger  added  mass  forces  mean  that  the  amplitude  of 
the  lift  oscillations  continue  to  increase,  though  the  mean  lift  stays  approximately  constant. 

To  analyze  the  time-varying  behavior  in  more  detail,  we  take  the  discrete  Fourier  transform  of 
the  lift  coefficient  signal  in  time  for  each  trial.  The  results  for  this  are  shown  in  Figure  4.10,  for 
pitching  amplitude  a  a  =  1°-  For  cxm  <  25,  we  observe  one  dominant  frequency  peak,  corresponding 
to  the  pitching  frequency  /*.  Even  though  the  undisturbed  wake  is  stable  at  a  =  25°,  we  still 
observe  some  frequency  content  near  the  almost-unstable  vortex  shedding  mode  across  all  pitching 
frequencies.  For  larger  angles  of  attack,  there  is  a  second  major  peak  in  the  spectra,  corresponding 
to  the  vortex  shedding  frequency  at  the  given  ocm •  In  the  region  where  these  two  frequencies 
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Figure  4.5:  Mean  drag  coefficient  for  a  range  of  pitching  amplitudes  ola,  frequencies,  /*,  and  mean 
angles  of  attack,  a«.  Pitching  is  about  the  midchord. 


are  similar,  there  appear  to  be  complex  interactions  between  the  dynamics  associated  with  each 
frequency.  As  olm  grows  larger,  a  distinct  peak  emerges  at  the  first  harmonic  of  the  fundamental 
vortex  shedding  frequency,  further  complicating  the  frequency  response  of  the  system,  which  can 
now  include,  at  very  least,  sums  and  differences  of  multiples  of  these  frequencies. 

4.3.3  Wake  analysis  with  dynamic  mode  decomposition 

Here  we  analyze  the  flow  field  in  the  wake  of  the  body,  in  an  attempt  to  study  the  underlying 
physics  behind  the  phenomena  observed  in  section  4.3.2.  We  will  make  use  of  the  dynamic  mode 
decomposition  (DMD)[125,  126]  which  is  a  technique  that  can  extract  dynamical  content  from  data, 
in  the  form  of  spatial  modes  and  their  associated  growth/decay  rates  and  frequencies  of  oscillation. 
We  refer  the  reader  to  recent  references  for  details  of  the  DMD  algorithm[126,  119,  149,  160]  and 
its  variants[27,  164,  71,  62,  63,  40]. 

We  begin  by  considering  a  base  angle  olm  =  20°,  which  is  prior  to  the  Hopf  bifurcation  at 
which  unforced  vortex  shedding  occurs,  and  thus  has  a  stable  equilibrium  wake.  Figure  4.11  shows 
vorticity  field  snapshots  for  a  variety  of  forcing  frequencies  with  amplitude  a  a  =  1°.  We  observe 
that  at  around  f*  ~  0.3,  the  vorticity  field  is  qualitatively  different,  with  periodic  vortex  shedding 
being  excited  by  the  pitching  motion.  This  immediately  suggests  that  the  increased  lift  observed 
at  these  pitching  frequencies  arises  due  to  enhanced  vortex  formation.  We  note  that  the  vorticity 
fields  for  pitching  about  the  leading  edge  (left)  and  midchord  (right)  show  similar  results,  with 
leading-edge  pitching  generally  leading  to  stronger,  more  distinct  vortices  forming  closer  to  the 
airfoil. 
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Figure  4.6:  Mean  of  lift-to-drag  ratio  for  a  range  of  pitching  amplitudes  a  a-,  frequencies,  /*,  and 
mean  angles  of  attack,  olm-  Pitching  is  about  the  midchord. 


The  vorticity  fields  shown  in  Figure  4.12  have  the  same  parameters  as  those  for  Figure  4.11, 
but  with  a  larger  pitching  amplitude  of  a  a  =  5°.  We  again  see  the  same  phenomena  where  the 
forcing  excites  vortex  shedding,  but  here  distinct  vortices  form  and  persist  downstream  over  a  wider 
range  of  frequencies.  This  is  consistent  with  the  results  from  section  4.3.2,  where  for  the  olm  =  20° 
case  in  Figures4.4  we  observe  an  enhanced  lift  over  a  wider  range  of  frequencies  for  higher  forcing 
amplitudes. 

To  investigate  further  the  dynamics  of  the  wakes,  we  perform  DMD  on  the  vorticity  field.  We 
use  191  snapshots  for  each  case,  with  a  time  gap  between  snapshots  of  At  =  0.1  c/U .  We  take  data 
from  a  region  downstream  of  the  body  as  our  domain  for  DMD,  to  avoid  complications  associated 
with  having  the  moving  body  in  the  domain.  For  clarity  of  results,  we  truncate  the  rank  of  our 
data  to  20  proper  orthogonal  decomposition  (POD)  [66]  modes.  This  is  found  to  capture  at  least 
99%  of  the  energy  of  the  flow. 

While  it  is  most  often  applied  to  data  from  systems  without  external  forcing,  DMD  can  be 
used  to  examine  forced  systems[147,  105,  92],  though  one  must  take  care  in  the  interpretation  of  its 
outputs.  In  this  context,  we  show  that  it  can  be  an  effective  tool  to  determine  both  the  frequency 
content  of  the  wake,  and  also  if  the  wake  “locks  on”  to  a  vortex  shedding  mode.  As  an  aside,  we 
note  that  more  rigorous  connections  between  DMD  and  Fourier  analysis  can  be  made[27]. 

Figure  4.13  shows  the  DMD  eigenvalues,  as  well  as  the  four  largest  amplitude  modes,  for  a  variety 
of  forcing  frequencies,  for  pitching  about  the  leading  edge  (left)  midchord  (right)  with  amplitude 
a  a  =  1°.  In  this  section,  we  restrict  our  attention  to  frequencies  near  the  vortex  shedding  frequency 
and  frequency  of  maximum  lift.  Results  at  lower  and  higher  frequencies  are  similar  to  the  minimum 
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Figure  4.7:  Mean  lift  coefficient  increment  over  the  fixed  airfoil  value,  normalized  by  the  amplitude 
of  pitching  a  a,  across  a  range  of  pitching  amplitudes  a  a,  frequencies,  /*,  and  mean  angles  of  attack, 
otM-  Pitching  is  about  the  midchord. 
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Figure  4.8:  Mean  lift  coefficient  increment  over  the  fixed  airfoil  value,  normalized  by  the  amplitude 
of  pitching  a  a-,  across  a  range  of  pitching  amplitudes  a  a,  frequencies,  /*,  and  mean  angles  of  attack, 
otM-  Pitching  is  about  the  leading  edge. 
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Figure  4.10:  Power  spectral  densities  of  the  lift  force  for  pitching  motion  about  the  midchord  with 
amplitude  oia  =  1°.  For  clarity,  the  spectra  for  each  forcing  frequency  is  shifted  by  one  order  of 
magnitude,  so  the  absolute  scale  of  the  vertical  axis  is  not  significant. 
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Figure  4.11:  Instantaneous  vorticity  fields  for  pitching  about  the  leading  edge  (left)  and  midchord 
(right)  at  a  variety  of  frequencies,  with  a.M  =  20°  and  a  a  =  1°. 
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Figure  4.12:  Instantaneous  vorticity  fields  for  pitching  about  the  leading  edge  (left)  and  midchord 
(right)  at  a  variety  of  frequencies,  with  a.M  =  20°  and  a  a  =  1°. 
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and  maximum  frequencies  considered  here.  Note  that  the  real  and  imaginary  components  of  the 
eigenvalues  (in  continuous  time)  indicate  the  growth  rate  and  frequency  of  the  associated  DMD 
modes. 

For  this  mean  angle  of  attack  (oim  =  20°),  the  maximum  increase  in  lift  occurrs  at  /*  =  0.3. 
At  this  frequency,  the  DMD  eigenvalues  are  at  the  frequency  of  pitching  and  its  harmonics.  Since 
vortex  shedding  is  clearly  observed  in  Figure  4.11,  this  suggests  that  it  is  locking  onto  the  forcing 
frequency.  For  the  surrounding  frequencies  (/*  =  0.25  and  0.35),  we  observe  that  the  eigenvalues 
are  scattered  across  a  range  of  frequencies.  This  suggests  that  vortex  shedding  is  not  fully  locking 
on  to  the  forcing  frequency,  and  that  there  is  a  broad  range  of  frequency  content  present.  We 
also  find  that  the  second,  third  and  fourth  DMD  modes  all  have  similar  spatial  wavelength,  again 
indicating  that  certain  features  persist  across  a  range  of  frequencies.  Note  that,  unlike  POD  modes, 
DMD  modes  are  not  required  to  be  orthogonal,  and  thus  there  can  exist  multiple  DMD  modes  that 
are  similar  to  each  other,  having  slightly  different  frequencies  of  oscillation. 

If  we  move  further  away  from  /*  =  0.3,  we  find  that  the  DMD  modes  again  lie  on  harmonics 
of  the  forcing  frequency.  In  these  cases  (i.e.,  for  /*  =  0.2  and  0.4  plotted  here,  and  for  other  cases 
not  shown),  it  appears  that  the  least  stable  mode  in  the  wake  is  not  being  excited,  as  we  are  not 
close  enough  in  frequency  to  do  so.  This  also  explains  why  there  is  no  significant  lift  increase  at 
these  frequencies.  The  similarity  between  pitching  at  the  leading  edge  and  midchord  suggests  that 
these  observations  are  not  particularly  sensitive  to  pitch-axis  location. 

Figure  4.14  shows  plots  of  the  same  quantities  and  parameters  as  Figure  4.13,  but  for  pitching 
with  amplitude  a  a  =  5°.  Broadly  speaking,  the  results  are  similar,  with  DMD  eigenvalues  lying  on 
the  imaginary  axis  at  multiples  of  the  forcing  frequency  for  f*  =  0.2,  0.3  and  0.4,  and  located  in 
a  significantly  more  scattered  arrangement  with  numerous  modes  of  similar  frequency  content  and 
spatial  structure  for  /*  =  0.25  and  0.35.  We  find  that  the  oscillating  modes  have  higher  amplitude 
in  this  case,  with  the  “lock-on”  mode  having  higher  amplitude  than  the  mean  mode  for  leading  edge 
pitching  at  /*  =  0.3.  We  further  note  that,  unlike  the  case  for  ola  =  1°,  the  mean  flow  (i.e.,  the 
DMD  mode  with  a  corresponding  eigenvalue  close  to  0)  looks  substantially  different  between  low 
frequency  oscillation,  where  the  wake  seems  to  spread  out  substantially  a  short  distance  downstream 
of  the  airfoil,  and  high  frequency  oscillation,  where  the  wake  appears  to  be  narrower  with  regions 
of  high  vorticity  persisting  far  downstream.  Referring  back  to  Figure  4.12,  the  enhanced  spreading 
of  the  mean  flow  is  seemingly  due  to  the  larger  vortices  drifting  above  (for  negative  vorticity)  and 
below  (positive  vorticity)  the  airfoil  as  they  are  shed.  This  phenomenon  could  also  explain  why 
the  lift-to-drag  ratio  is  often  maximized  at  pitching  frequencies  slightly  larger  than  those  at  which 
maximum  lift  occurs.  The  enhanced  spreading  at  lower  frequencies  increases  drag,  so  lift-to-drag 
benefits  exist  only  at  higher  frequencies,  where  the  wake  remains  thinner. 

Thus  far,  we  have  only  analyzed  the  wake  for  a  case  where  the  base  flow  is  stable.  We  now  turn 
our  attention  to  the  case,  olm  =  30°,  where  the  flow  exhibits  a  vortex  shedding  limit  cycle.  Figure 
4.15  shows  vorticity  field  snapshots,  as  well  as  DMD  eigenvalues  and  modes  (only  showing  the  four 
highest  amplitude  modes)  for  pitching  about  a  mean  angle  au  =  30°  with  amplitude  a  a  =  1°,  for  a 
variety  of  pitching  frequencies.  The  wakes  for  f*  =  0.2,  0.3,  0.35  and  0.4  all  look  very  similar,  with 
high  amplitude  DMD  eigenvalues  at  the  natural  vortex  shedding  frequency,  and  lower  amplitude 
eigenvalues  at  the  pitching  frequency.  This  is  consistent  with  Figure  4.10,  where  the  subplot  for 
qm  =  30°  shows  a  larger  peak  at  the  vortex  shedding  frequency,  and  a  lower  peak  at  the  forcing 
frequency.  For  f*  =  0.25,  the  forcing  and  vortex  shedding  frequencies  are  almost  identical,  which 
appears  to  give  stronger  vortex  formation,  and  in  turn  the  increase  in  lift  seen  in  Figure  4.4. 
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Figure  4.13:  DMD  eigenvalues,  and  the  real  components  of  the  4  largest  amplitude  modes  (ordered 
by  amplitude),  for  pitching  about  the  leading  edge  (left),  and  midchord  (right)  at  a  variety  of 
frequencies,  with  amplitude  a  a  =  1°,  and  mean  angle  of  attack  au  =  20°.  Eigenvalues  that  are 
colored  red  correspond  to  modes  with  larger  amplitudes. 
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Figure  4.14:  DMD  eigenvalues,  and  the  real  components  of  the  4  largest  amplitude  modes  (ordered 
by  amplitude),  for  pitching  about  the  leading  edge  (left),  and  midchord  (right)  at  a  variety  of 
frequencies,  with  amplitude  a  a  =  5°,  and  mean  angle  of  attack  au  =  20°.  Eigenvalues  that  are 
colored  red  correspond  to  modes  with  larger  amplitudes. 
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Figure  4.15:  Wake  snapshots,  DMD  eigenvalues,  and  the  real  components  of  the  4  largest  amplitude 
modes  (ordered  by  amplitude),  for  midchord  pitching  at  a  variety  of  frequencies,  with  amplitude 
a  a  =  1°,  and  mean  angle  of  attack  cym  =  30°.  Eigenvalues  that  are  colored  red  correspond  to 
modes  with  larger  amplitudes. 
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4.4  Discussion  and  conclusions 

For  angles  of  attack  below  the  critical  angle  where  unforced  vortex  shedding  first  occurs  (ac) ,  it  was 
found  that  pitching  at  a  certain  frequency  can  excite  vortex  shedding  in  the  wake,  leading  to  higher 
mean  lift.  The  magnitude  and  width  in  frequency  range  of  the  lift  increment  increases  significantly 
as  the  forcing  amplitude  increases.  This  is  perhaps  a  similar  phenomenon  to  the  widening  of  the 
“resonance  horn”  [20]  observed  by  Choi  et  al.  [29]  for  the  case  of  a  surging  and  plunging  airfoil 
(though  in  that  case,  the  system  was  above  the  critical  bifurcation  parameter).  Note,  however, 
that  comparing  Figures  4.4  and  4.14  indicates  that  there  can  be  a  significant  increase  in  lift  even 
without  frequency  lock-on. 

For  oiM  above  the  critical  angle  for  vortex  shedding,  there  is  a  similar  peak  in  the  mean  lift 
coefficient  when  the  pitching  frequency  is  close  to  the  natural  vortex  shedding  frequency  or  its  first 
harmonic.  When  the  natural  and  forcing  frequencies  are  different,  the  interactions  between  the  two 
frequencies  can  lead  to  complex  frequency  spectra  in  the  forces  and  wakes.  Note  in  particular  that 
pitching  some  amount  below  the  vortex  shedding  frequency  can  lead  to  a  notable  decrease  in  mean 
lift  for  a  a  >  5°  at  qm  >  35°. 

While  periodic  pitching  at  the  preferred  frequency  where  vortex  shedding  is  excited  or  enhanced 
also  leads  to  an  increase  in  drag,  the  differences  in  effect  that  pitching  has  on  the  two  force  compo¬ 
nents  leads  to  an  increase  in  the  lift-to-drag  ratio  for  frequencies  slightly  above  the  frequency  for 
which  the  lift  is  maximized.  It  is  interesting  to  note  that  while  the  frequency  at  which  the  maximum 
lift  increment  occurs  decreases  slightly  as  increases  (in  agreement  with  the  slight  decrease  in 
natural  vortex  shedding  frequency) ,  the  frequency  giving  maximum  lift-to-drag  ratio  increases  with 
aM-  Indeed,  for  high  q,m  it  seems  that  the  maximum  lift-to-drag  ratio  typically  occurs  between 
the  peaks  in  lift  and  drag  located  at  the  vortex  shedding  frequency  and  its  first  harmonic,  where 
there  is  a  local  lift  minimum,  and  thus  also  a  slightly  more  substantial  drag  minimum. 

Beyond  ac,  the  interactions  between  the  pitching  and  vortex  shedding  frequencies  can  lead  to 
complex  frequency  spectra  in  both  the  forces  (as  seen  in  Figure  4.10)  and  wakes  (Figures  4.13- 
4.15).  There  are  numerous  methods  by  which  one  can  analyze  frequency  content,  and  here  we  show 
how  DMD  can  clearly  distinguish  between  cases  where  all  frequencies  present  are  harmonics  of  the 
pitching  frequency,  and  where  there  is  a  broad  range  of  frequency  content. 

There  has  been  much  effort  in  the  past  to  understand,  model,  and  predict  unsteady  aerody¬ 
namic  forces,  moments,  pressures,  and  indeed  many  other  quantities  of  interest  for  moving  airfoils. 
Particularly  in  the  case  of  separated  flow,  these  dynamics  are  often  highly  nonlinear.  One  might 
seek  to  get  around  this  by  linearizing  about  a  certain  fixed  point  (say  an  angle  of  attack),  in  the 
hope  that  a  linear  model  would  at  least  be  locally  accurate.  The  findings  presented  here  suggest 
that  such  an  approach  might  be  problematic,  since  a  linear  model  (e.g.,  a  transfer  function)  can 
only  predict  the  magnitude  and  phase  of  a  response  to  sinusoidal  forcing,  but  not  any  change  in 
the  mean  value.  Thus,  such  effects  must  be  accounted  for  separately,  or  a  more  complex  modeling 
framework  used. 

Further  work  in  this  investigation  will  compare  the  findings  to  the  results  of  stability  analyses 
of  the  stable  and  unstable  equilibrium  wakes,  and  the  mean  of  the  vortex  shedding  wake,  and  will 
use  wind  tunnel  experiments  to  investigate  whether  the  phenomena  observed  persist  for  higher 
Reynolds  numbers. 
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Chapter  5 

Nonlinear  reduced-order  models  of 
fluids  systems  using  extended 
dynamic  mode  decomposition 


Data-driven  approximations  to  the  Koopman  operator  have  much  potential  for  capturing  and  il¬ 
luminating  the  dynamics  exhibited  by  the  Navier-Stokes  equations.  In  this  work,  we  show  that 
the  elements  of  an  identified  finite  dimensional  approximation  to  the  Koopman  operator  can  be 
utilized  for  the  construction  of  accurate  nonlinear  reduced-order  models.  We  present  a  modification 
to  the  extended  dynamic  mode  decomposition  algorithm  through  the  inclusion  of  a  regularization 
parameter,  which  we  find  often  gives  more  accurate  models.  The  performance  of  models  identified 
using  our  proposed  method  are  compared  to  those  found  by  performing  a  Galerkin  projection  of 
the  governing  equations  onto  proper  orthogonal  decomposition  modes,  for  the  canonical  case  of 
two-dimensional  flow  past  a  circular  cylinder.  We  demonstrate  that  identifying  nonlinear  models 
using  EDMD  is  particularly  advantageous  when  the  data  available  is  noisy,  or  is  only  available 
within  certain  regions  of  space  or  time. 

5.1  Introduction 

A  much-desired  goal  in  fluid  mechanics,  and  indeed  many  other  fields,  is  to  obtain  simple  models 
that  are  capable  of  predicting  the  behavior  of  seemingly  complex  systems.  Low-dimensional  models 
can  not  only  improve  our  fundamental  understanding  of  such  systems,  but  are  often  required  for 
purpose  of  efficient  and  accurate  prediction,  estimation  and  control.  Broadly  speaking,  one  can 
obtain  low-dimensional  information  about  a  system  (whether  it  be  in  the  form  of  a  reduced-order 
model,  or  simply  spatial  modes  corresponding  to  certain  energetic  or  dynamic  characteristics) 
in  numerous  ways,  potentially  using  some  combination  of  data  collected  from  simulations  and 
experiments,  and  theoretical  knowledge  of  the  system,  such  as  the  governing  partial  differential 
equations  (PDEs). 

Purely  data-driven  methods  can  include  those  developed  particularly  for  fluids  applications, 
such  as  the  dynamic  mode  decomposition  (DMD)  [125,  126],  or  those  which  are  appropriated  from 
other  communities,  such  as  the  eigensystem  realization  algorithm  (ERA)  [82,  72],  which  was  first 
applied  to  study  spacecraft  structures,  but  has  more  recently  been  appropriated  to  model  a  wide 
range  of  fluids  systems  [25,  2,  68,  69,  22,  23,  15,  67,  49].  DMD  has  been  used  to  study  a  wide  range 
of  problems  arising  in  fluid  mechanics  (e.g.,  [119,  128,  127,  100,  129,  43,  58,  70,  83,  93,  55,  123,  45]), 
with  many  subsequent  variants  further  increasing  its  range  of  utility  and  applications  [27,  164,  149, 
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71,  62,  63,  40,  106,  124,  98,  122],  One  limitation  of  the  methods  mentioned  so  far  is  that  (when 
considered  in  the  context  of  data-driven  reduced-order  modeling  techniques)  they  are  linear,  in  the 
sense  that  the  reduced  order  model  that  is  identified  is  in  the  form  of  a  linear  system  of  ordinary 
differential  equations  (ODEs). 

While  there  have  been  a  number  of  examples  of  nonlinear  data-driven  modeling  techniques 
used  in  fluids  applications  [85,  102,  104,  52,  11,  37,  77,  60,  39,  24],  their  widespread  use  has  been 
more  limited,  and  the  underlying  theory  is  less  established,  than  linear  techniques.  More  details 
concerning  the  application  of  data-driven  modeling  techniques  in  fluid  mechanics  can  be  found  in 
recent  review  articles  [21,  113]. 

Perhaps  the  most  common  method  to  obtain  a  nonlinear  reduced  order  model  for  fluids  systems 
comes  via  a  projection  of  the  governing  equations  onto  a  low-dimensional  basis  that  is  optimal  for 
capturing  the  energy  of  the  data,  i.e.,  the  proper  orthogonal  decomposition  (POD)  [86,  17,  66],  a 
procedure  referred  to  as  Galerkin  projection.  Galerkin  projection  (GP)  has  been  used  to  extract 
models  for  many  different  fluids  systems,  a  non-exhaustive  list  includes  flow  past  a  cylinder  at 
low  Reynolds  number  [41,  96,  95],  grooved  channels  [41]  the  wall  region  of  turbulent  boundary 
layers  [6,  103],  flat  plate  boundary  layers  [111],  turbulent  plane  Couette  flow  [91,  134],  turbulent 
pipe  flow  [19]  cavity  oscillations  [116,  115],  mixing  layers  [110,  150,  13],  and  compressible  flows 
[118].  One  significant  drawback  of  GP  models  is  that  they  ignore  modes  that  are  low  in  energy, 
but  are  required  for  the  dissipation  of  energy  in  the  full  system.  A  number  of  modifications  have 
been  proposed  to  address  this  concern,  as  well  as  other  issues  with  such  models.  Refs.  [6]  and 
[103]  use  an  eddy  viscosity  term  that  accounts  for  energy  dissipated  into  unmodeled  modes,  [99] 
investigates  a  hierachy  of  eddy  viscosity  formulations,  while  [157,  158]  incorporate  LES  closure 
modeling  strategies.  Refs.  [35]  and  [36]  summarize  a  number  of  calibration  techniques  that  can 
be  used  to  improve  the  accuracy  of  Galerkin  models,  and  also  discuss  the  various  ways  in  which 
the  error  of  such  models  can  be  quantified.  Ref.  [12]  employs  a  subspace  rotation  technique  to 
stabilize  the  models,  which,  unlike  other  calibration  techniques,  maintains  consistency  with  the 
original  governing  equations.  Ref.  [13]  imposes  constraints  to  balance  the  turbulent  kinetic  energy 
of  the  resulting  model.  All  of  these  modifications  of  Galerkin  projection  increase  the  “data-driven” 
nature  of  the  method.  Ref.  [97]  gives  an  in-depth  summary  and  analysis  of  many  issues,  variations, 
progress,  and  open  problems  on  the  topic  of  Galerkin  projection  models,  while  [84]  gives  a  clear 
expository  introduction  of  the  main  ideas  in  Galerkin  projection,  with  examples. 

While  we  mentioned  above  that  DMD  could  be  classified  as  a  “linear”  method,  connections 
between  the  DMD  algorithm  and  the  Koopman  operator  [119,  89]  give  promise  that  it  can  ultimately 
be  used  to  model  and  understand  nonlinearities.  In  this  work,  we  will  propose  a  method  for 
obtaining  nonlinear  reduced  order  models  that  is  based  upon  a  recently-developed  extension  of 
DMD  [160],  referred  to  as  extended  DMD  (EDMD),  in  which  nonlinearities  can  be  accounted  for 
by  an  appropriate  choice  of  observables.  In  particular,  we  will  explore  how  this  algorithm  can  be 
tailored  to  identify  nonlinear  models  that  have  the  same  (or  similar)  form  as  those  that  would  be 
obtained  through  projection  of  the  governing  equations  onto  a  low-dimensional  basis  obtained  from 
data.  In  this  sense,  the  nonlinear  models  that  we  identify  in  this  work  will  come  from  EDMD,  but 
will  have  similar  form  to  those  given  by  GP.  We  note  that  this  is  a  different  approach  to  using 
dominant  DMD  modes  as  an  alternative  to  POD  modes  for  a  basis  for  projection  [143].  A  further 
recently-proposed  approach  uses  DMD  to  efficiently  approximate  just  the  nonlinear  component  of 
the  dynamics  [5],  as  an  alternative  to  the  discrete  empirical  interpolation  method  [14,  26]. 

We  lastly  note  that  several  previous  works  consider  this  dichotomy  between  using  projection 
onto  governing  equations,  or  using  time-resolved  data,  to  identify  models,  most  often  for  linear 
systems.  For  example,  [69]  considers  models  identified  directly  using  ERA  and  models  identified 
by  considering  each  of  the  pertinent  physical  processes  individually,  while  [64]  discusses  the  use  of 
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both  GP  models  and  ARMAX  system  identification  methods  with  physically  motivated  terms  to 
model  (and  subsequently  control)  flow  over  a  backwards-facing  step. 

The  structure  of  this  work  is  as  follows.  Sect.  5.2  introduces  POD  and  Galerkin  projection. 
Aside  from  GP  being  used  for  the  purposes  of  comparison  to  EDMD  models,  the  POD  section  will 
also  be  important  for  the  purposes  of  defining  observables  for  EDMD.  Sect.  5.3  describes  DMD  and 
EDMD,  and  introduces  the  particular  form  of  the  EDMD,  as  well  as  a  regularized  variant  thereof, 
that  will  be  utilized  for  system  identification  purposes.  Models  will  be  identified  and  tested  in 
Sect.  5.4  on  data  obtained  from  numerical  simulations  of  flow  past  a  circular  cylinder,  the  results 
of  which  are  discussed  in  further  detail  in  Sect.  5.5. 


5.2  Proper  orthogonal  decomposition  and  Galerkin  projection 

While  the  main  focus  of  this  work  concerns  reduced-order  models  obtained  from  EDMD  algorithms, 
GP  models  will  serve  as  both  a  basis  for  comparison  of  model  performance,  and  to  guide  our  choice 
of  observable  functions  when  using  EDMD.  With  this  in  mind,  in  this  summary  we  give  a  brief 
summary  of  POD  and  Galerkin  projection. 


5.2.1  Proper  orthogonal  decomposition 

The  goal  of  the  proper  orthogonal  decomposition  (POD)  is  to  obtain  a  set  of  empirical  spatial 
modes  that  optimally  represent  a  given  dataset  from  an  energetic  standpoint.  Assume  that  we  can 
decompose  the  dynamics  of  some  system  u(x,t )  (which  could  be  the  time- varying  velocity  field  of 
a  fluid,  say)  by 

OO 

u(x,  t )  =  U0(x)  +  ^2  ui(x)ai(t),  (5.2.1) 

i=i 

where  uq(x)  is  some  fixed  (often  average)  data,  and  {■Uj(x)}?fi1  are  a  set  of  orthonormal  basis 
functions  (modes).  POD  takes  these  modes  to  be  those  which  successively  capture  the  most  energy 
of  the  velocity  field.  Each  POD  mode  ti;  satisfies  the  integral 


Rij(x,  x')ui(x')dx'  =  A Ui(x), 


where  R,hj(x,  y)  =  E[ui(x)uj(y)\,  with  E  being  the  expectation.  As  indicated  by  Eq.  (5.2.1), 
POD  is  normally  performed  after  first  subtracting  the  mean  (or  perhaps  an  equilibrium  point) 
from  the  data.  This  approach  has  the  advantage  that  uq  satisfies  the  required  non-homogeneous 
boundary  conditions,  meaning  that  all  other  modes  Ui  will  satisfy  homogenous  boundary  conditions, 
so  any  linear  combination  of  modes  of  the  form  given  by  Equation  (5.2.1)  will  automatically  satisfy 
the  correct  boundary  conditions  of  the  problem  at  hand.  In  discrete  terms,  if  we  arrange  finite¬ 
dimensional  data  collected  from  a  simulation  or  experiment  into  a  matrix  Y,  with  each  column 
representing  a  snapshot  of  data  at  a  given  time,  then  the  POD  modes  are  the  columns  of  U  in  the 
singular  value  decomposition  Y  =  UYJV* .  Here  the  ith  entry  of  the  diagonal  matrix  X  corresponds 
to  the  energy  contained  in  the  ith  POD  mode.  In  this  discrete  formulation,  for  simplicity  we  are 
omitting  any  rescaling  of  the  data  that  should  be  performed  so  that  the  modes  are  orthonormal 
with  respect  to  the  usual  inner  product.  That  is,  if  u%  and  uj  are  columns  of  U,  then  we  really 
should  have 

n 

/  u*(x')u.i(x')dx'  «  y^  u*(xk)ui(xk)dxk  =  Sij, 
k=  1 
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rather  than  u*Uj  =  Sjj.  The  original  data  Y  can  then  be  represented  in  terms  of  POD  coefficients 
by  Y  =  U*Y .  If  we  wish  to  reduce  the  dimension  of  this  data,  we  may  do  so  in  an  optimal  way 
(with  respect  to  energy  content)  by  simply  truncating  the  columns  of  U  beyond  a  certain  point, 
which  corresponds  to  removing  POD  modes  that  are  of  sufficiently  low  energy.  Doing  this  gives 
a  reduced  order  approximation  of  the  data  Yr  =  U*Y,  where  Ur  contains  the  first  r  columns  of 
U.  Note  that  there  are  alternative  truncation  techniques  that  may  be  more  effective  than  energy 
maximization  for  certain  applications,  for  example  balanced  POD  [112]  gives  a  reduced  order  linear 
state  space  model  that  is  optimal  with  respect  to  a  given  set  of  sensors  and  actuators. 

5.2.2  Galerkin  projection 

The  idea  behind  GP  is  to  approximate  the  governing  PDEs  that  describe  a  given  system  with  a 
low- dimensional  set  of  ODEs.  This  is  accomplished  by  projecting  the  equations  onto  spatial  POD 
modes  identified  using  the  methods  described  in  Sect.  5.2.1.  We  begin  with  the  incompressible 
Navier-Stokes  equations: 

Qu 

— -  =  —(u  ■  V)u  —  uAu  —  VP 

dt  y  J  (5.2.2) 

V  •  u  =  0. 

If  we  take  the  (spatial)  inner  product  of  Eq.  (5.2.2)  with  a  given  mode  Uj,  we  obtain 

=  —  (( u  ■  V)u,  Uj)  —  v  (A u,  Uj)  —  (VP,  Uj) .  (5.2.3) 

Substituting  in  the  finite-dimensional  approximation  of  Eq.  (5.2.1),  we  obtain 

a  =  La  +  Q(a,  a)  +  f,  (5.2.4) 

were  L  is  a  linear  operation  (i.e.,  a  matrix),  Q  is  a  bilinear  operator  (which  can  be  represented  as 
a  3-tensor),  and  f  is  a  vector,  each  defined  based  on  the  identified  spatial  POD  modes  by 

Lij  =  -v  (Auj,Ui) ,  Qijk  =  ~  (( Uj  ■  V)uk,Ui) ,  f  =  -  {Vp,  Uj)  .  (5.2.5) 

This  gives  a  means  of  approximating  the  Navier-Stokes  equations  by  a  set  of  nonlinear  ODEs.  As 
mentioned  in  Sect.  6.1,  there  are  many  modifications  that  have  been  proposed  for  this  general 
procedure,  most  typically  to  account  for  the  energy  transfer  to  unmodeled  modes  (i.e.,  the  energy 
cascade  to  finer  spatial  scales).  For  cases  where  spatial  symmetries  exist  (e.g.,  in  the  streamwise 
and  azimuthal  directions  for  circular  pipe  flow),  one  can  show  that  the  POD  modes  must  become 
Fourier  modes,  which  can  simplify  their  computation.  It  is  also  possible  to  “factor  out”  such 
symmetries  by  using  an  optimally  chosen  moving  frame  of  reference [114,  117]. 


5.3  Reduced-order  models  using  extended  dynamic  mode  decom¬ 
position 

This  section  introduces  our  proposed  modeling  approach.  Sect.  6.2.1  describes  the  DMD  algorithm, 
before  Sect.  5.3.2  discusses  the  EDMD  extension  and  how  it  may  be  used  to  obtain  nonlinear 
models.  Sect.  5.3.3  gives  a  regularized  modification  to  EDMD  that  we  find  to  be  advantageous 
when  using  EDMD  for  such  purposes.  When  viewed  as  a  method  for  reduced-order  modeling,  the 
main  difference  between  the  approach  discussed  here  and  that  introduced  in  Sect.  5.2  is  in  how  the 
temporal  dynamics  are  identified:  GP  uses  the  governing  equations,  whereas  DMD/EDMD  uses 
only  data  to  identify  dynamics. 
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5.3.1  Dynamic  Mode  Decomposition 

u 

Following  the  notation  introduced  in  [113],  suppose  we  collect  pairs  of  snapshots  of  data  (y,;,  y,-  ), 
which  are  separated  by  a  fixed  time  interval  At.  y,  is  often  taken  to  be  a  velocity  field  snapshot,  but 
more  generally  it  is  a  vector  consisting  of  a  number  of  observations/measurements  of  an  evolving 
dynamical  system  at  a  given  point  in  time.  We  arrange  the  snapshot  pairs  as  columns  of  the 
matrices  Y  and  Y#  (where  the  data  in  a  given  column  of  Y*  was  collected  At  after  the  equivalent 
column  of  Y).  Note  that  if  we  begin  with  a  sequential  time-series  of  data,  then  Y  and  Y #  are 
simply  formed  using  all  of  the  snapshots  excluding  the  last  and  first  columns,  respectively.  Let  Y 
and  Y#  each  be  n  by  m  matrices,  so  we  have  m  pairs  of  snapshots,  each  of  size  n.  One  can  define 
(see  [149])  DMD  modes  and  eigenvalues  as  the  eigendecomposition  of  the  matrix 

A  =  Y*Y+,  (5.3.1) 


where  Y+  is  the  Moore-Penrose  pseudoinverse  of  a  matrix  Y .  While  useful  as  a  definition,  one  does 
not  typically  compute  A  directly,  since  in  the  typical  case  where  n  >  m,  4  is  a  very  large  n  x  n 
matrix  with  rank  at  most  min(m,  n).  Thus  it  is  usually  beneficial  to  first  project  the  data  onto 
a  lower  dimensional  subspace,  which  is  often  taken  from  the  columns  of  U  in  the  singular  value 
decomposition  Y  =  LTSV*  (possibly  truncating  columns  of  U  corresponding  to  small  singular 
values).  From  Sect.  5.2.1,  this  subspace  consists  of  the  POD  modes  of  the  data.  If  we  let  y  =  U*yk 
be  the  representation  of  a  given  snapshot  y  in  the  POD  basis,  and  let  Y  =  U*Y  and  Y#  =  U*Y # 
following  the  same  notation  as  introduced  in  Sect.  5.2.1,  then  the  equivalent  of  Eq.  (6.2.1)  in  POD 
space  is 

A  =  Y*Y+.  (5.3.2) 


A  reduced  order  approximation  of  the  identified  dynamics  can  be  obtained  by  restricting  the  number 
of  POD  modes  upon  which  to  project.  This  is  achieved  by  taking  the  first  r  columns  of  U,  and 
working  in  the  space  of  POD  coefficients  =  U* y^.  DMD  then  results  in  a  simple  linear  dynamical 
system  model  of  the  form 

ak+i  —  (5.3.3) 


where  Ar  =  U*rY*Y+Ur. 


5.3.2  Extended  dynamic  mode  decomposition  and  nonlinear  models 

A  fundamental  limitation  to  models  extracted  from  DMD  is  their  linearity,  which  can  make  them 
entirely  unable  to  model  fundamentally  nonlinear  features  that  arise  in  fluid  flows.  In  the  context 
of  DMD  being  an  approximation  of  the  Koopman  operator,  this  limitation  amounts  to  Koopman 
eigenfunctions  not  lying  within  the  span  of  the  data.  Following  [160],  rather  than  applying  the 
DMD  algorithm  directly  to  the  state  y,  we  may  define  to  define  a  set  of  observables  tp( y)  (where 
'0:  Mn  — >  Mk  for  some  k)  that  span  a  space  more  conducive  to  approximating  the  true  dynamics. 
Explicitly,  we  may  proceed  with  the  same  DMD  procedure  described  in  Sect.  6.2.1,  but  take  the 
columns  of  the  matrices  Y  and  Y #  to  be  pairs  (V'(yi),  V’(yj^))-  The  EDMD  matrix  (which  can  be 
viewed  as  a  finite-dimensional  approximation  to  the  Koopman  operator)  is  then  given  by 


A  =  Y*Y+  = 


V’(yf)  V>(yf )  •  •  •  V’(ym)  W>(yi)  V>(y2)  ■  ■  ■  V>(ym)r 


7#^ 

'2 


(5.3.4) 


While  '0(Y)  =  U*Y  gives  an  optimal  transformation  of  the  data  from  an  energetic  perspective  (and 
is  what  Eq.  (6.2.2)  represents  in  the  EDMD  framework),  it  might  not  be  a  suitable  transformation 
for  correct  identification  of  the  dynamics.  To  this  end,  we  use  the  form  of  Galerkin  projection  models 
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to  guide  an  alternative  transformation.  Expressing  the  POD  coefficients  as  before  by  a 
let 


q  =  V’(y) 


a 

vec  (a  <8>  a) 


U*y,  we 
(5.3.5) 


Here  vec  (a  (g>  a)  denotes  a  vector  of  all  non-redundant  quadratic  couplings  between  POD  coeffi¬ 
cients,  i.e.  (ai)2,  aia2,  (c 12 )2,  etc.  If  we  keep  r  POD  modes,  then  we  have  observables  of  dimension 
k  =  r  +  r(r  +  l)/2,  with  the  possible  addition  of  an  additional  constant  observable  to  account  for 
the  mean  of  the  data.  We  will  not  closely  concern  ourselves  with  how  closely  we  may  approximate 
the  true  Koopman  operator  using  such  a  choice  of  observables,  but  will  rather  show  that,  in  any 
case,  the  elements  of  the  identified  dynamics  on  this  space  of  observables  may  be  used  to  construct 
a  nonlinear  model  of  the  system  dynamics.  To  this  end,  we  start  by  explicitly  writing  the  identified 
dynamics  using  this  approach  by 
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(5.3.6) 


Focusing  on  a  single  POD  coefficient  (i.e.,  one  of  the  first  r  rows  of  Eq.  5.3.6),  we  have  that 


ak{t  +  At)  =  °i(*)ai(*)>  (5-3'7) 

i=l  =  l 

j<i 


which  is  the  same  form  as  Eq.  (5.2.5),  but  in  discrete-time.  Thus,  by  taking  the  first  r  rows  of  the 
A  matrix  obtained  from  performing  EDMD  with  observables  given  by  Eq.  (5.3.5),  we  may  obtain 
a  system  of  nonlinear  equations  that  can  accurately  model  the  evolution  of  the  POD  coefficients  of 
the  system  in  discrete  time. 

While  we  may  expect  that  the  first  r  rows  of  A  give  an  accurate  model  for  the  evolution  of  POD 
coefficients,  in  general  we  should  not  expect  the  same  for  the  evolution  of  the  quadratic  monomials 
of  POD  coefficients.  This  is  because  the  equations  for  the  evolution  of  these  terms  should  involve 
cubic  terms,  which  are  not  spanned  by  our  observables.  Incidentally,  this  suggests  that  using  any 
basis  of  polynomial  observables  for  approximating  the  Koopman  operator  for  the  Navier-Stokes 
problems  might  be  problematic.  One  possible  alternative,  not  explored  in  this  work,  is  to  use  the 
Kernel  variant  of  EDMD  [161]  with  an  appropriately  chosen  kernel  function  that  better  spans  the 
Koopman  eigenfuctions. 


5.3.3  A  modification  to  DMD/EDMD 

When  in  the  least  squares  regime  (i.e.,  when  there  are  fewer  observables  than  the  number  of 
snapshot  pairs),  Eqs.  (6.2.1)  and  (5.3.4)  can  give  solutions  with  large  entries  in  A.  Empirically, 
these  entries  can  be  significantly  larger  in  magnitude  than  those  expected,  say,  from  performing 
a  Galerkin  projection.  It  also  appears  that  this  “overfitting”  can  give  models  that  lack  stability. 
To  mitigate  these  observations,  we  propose  a  simple  modification  to  DMD/EDMD  that  penalizes 
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the  size  of  entries  in  A.  Note  that  the  same  process  can  be  applied  to  both  DMD  and  EDMD, 
whenever  the  data  matrix  Y  has  more  columns  than  rows.  When  this  is  the  case,  Eq.  (5.3.4)  gives 
the  solution  to  the  minimization  problem 

A  =  &rgmmM\\MY  —  Y#\\2f.  (5.3.8) 

We  may  add  a  penalization  on  the  size  of  the  entries  of  A  to  formulate  a  joint  least  squares  problem 

A  =  argminM  (j| MY  -  Y*\\2F  +  /32||M||^  ,  (5.3.9) 

where  /3  is  a  parameter  that  determines  the  extent  to  which  large  entries  in  M  are  penalized.  The 
element-wise  nature  of  the  Frobenius  norm  means  that  Eq.  (5.3.9)  may  be  rearranged  to  give 

A  =  argminM  M  [Y  /3I]  -  \y*  oj  2  ,  (5.3.10) 

L  J 

where  I  and  0  are  appropriately  sized  matrices.  Note  that  this  modification  of  DMD/EDMD  is 
similar  to  the  principle  of  Tikhonov  regularization.  Eq.  (5.3.10)  has  an  explicit  solution,  given  by 

A  =  Y*[Y/3I}+.  (5.3.11) 

We  remark  that  this  regularization  is  equivalent  as  adding  to  the  data  set  pairs  of  snapshots  where 
each  observable  goes  from  some  nonzero  value  (3  to  0  over  a  timestep  of  At.  It  should  also  be 
noted  that  Tikhonov  regularization  methods  have  been  used  previously  as  a  method  of  calibrating 
GP  models [35],  and  have  also  been  used  for  other  system  identification  techniques,  such  as  finite 
impulse  response  models  [76].  More  generally,  one  can  add  a  variety  of  different  forms  of  penalty 
terms  to  obtain  a  desired  balance  between  competing  objectives.  For  example,  in  the  context  of 
EDMD,  [162]  uses  a  L\ $  minimization  to  obtain  what  in  our  notation  would  be  a  solution  with  a 
small  number  of  nonzero  columns  of  A. 

5.4  Example:  flow  past  a  circular  cylinder 

We  test  our  proposed  method  using  the  much-studied  example  of  2D  flow  past  a  circular  cylinder. 
Beyond  a  critical  Reynolds  number  of  approximately  47  [107],  the  equilibrium  becomes  unstable 
and  the  system  will  instead  converge  upon  a  limit  cycle  characterized  by  periodic  vortex  shedding. 
We  take  data  for  Re  =  60,  with  the  initial  condition  close  to  the  unstable  equilibrium.  The  data 
captures  the  initial  growth  of  an  instability  near  the  unstable  equilibrium,  through  to  convergence 
to  the  limit  cycle.  This  is  an  example  where  regular  DMD  will  fail  (in  the  sense  of  identifying  an 
appropriate  reduced  order  model),  since  the  process  (which  is  a  Hopf  bifurcation)  is  fundamentally 
nonlinear.  In  particular,  beyond  the  critical  Reynolds  number,  the  nonlinear  terms  must  become 
non-negligible  to  balance  the  growth  of  the  unstable  linear  dynamics,  leading  to  the  observed  limit 
cycle  behavior. 

We  will  explore  the  performance  of  both  EDMD  and  GP  models  on  this  system  for  a  range  of 
data,  including  that  which  is  noisy  and  spatially  truncated  or  sparse.  The  data  was  obtained  from 
direct  numerical  simulation  using  an  immersed  boundary  projection  method  [139,  34],  Selective 
frequency  damping  [3]  was  used  to  obtain  the  unstable  equilibrium  solution.  To  focus  on  the 
transitional  region  of  the  dynamics,  the  snapshots  to  be  used  were  collected  after  first  running  the 
simulation  from  the  unstable  equilibrium  for  250  jj  time  units.  The  simulation  was  performed  on  a 
domain  consisting  of  five  nested  grids,  as  shown  in  Fig.  5.1.  Each  grid  is  uniform,  with  the  finest 
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Figure  5.1:  Computational  domain  used  for  numerical  simulation  of  flow  past  a  circular  cylinder. 
Dashed  lines  represent  the  borders  of  each  nested  grid.  The  gray  shading  denotes  the  region  from 
which  data  was  collected  for  modeling  purposes 


grid  having  a  grid  spacing  of  0.02D  in  each  direction,  and  each  successively  larger  grid  having 
double  the  grid  spacing  of  the  previous.  The  full  domain  has  size  256 D  x  64 D.  This  large  domain 
was  chosen  so  as  to  resolve  both  the  flowfield  on  the  region  of  the  grid  used  for  analysis  (shown  in 
gray),  and  the  forces  incident  on  the  cylinder,  to  a  high  degree  of  accuracy. 

This  comparison  between  the  performance  of  EDMD  and  GP  models  begins  in  a  scenario  where 
both  models  work  relatively  well:  where  clean  data  encompassing  a  large  spatial  domain  is  available 
across  a  window  of  time  spanning  all  of  the  distinct  dynamic  regimes.  While  we  find  that  GP  models 
can  indeed  outperform  the  EDMD  models  that  we  identify  in  favorable  conditions,  we  will  proceed 
to  demonstrate  that  in  the  cases  where  the  data  is  noisy  or  restricted,  EDMD  models  can  perform 
substantially  better. 

5.4.1  Data  arrangement  and  selection 

The  data  to  be  used  was  taken  from  the  gray  shaded  region  in  Fig.  5.1.  As  mentioned  previously, 
we  evolve  the  system  from  the  equilibrium  for  some  time  (250  jj)  before  collecting  data,  then  collect 
1000  snapshots  separated  by  a  uniform  timestep  At  =  0.2^.  This  time  interval  spans  the  growth 
of  the  instability  from  near  the  equilibrium,  through  to  the  convergence  of  the  flow  to  a  periodic 
limit  cycle.  We  note  that  Galerkin  projection  in  particular  is  quite  sensitive  to  the  resolution  and 
extent  of  data  chosen.  For  this  reason,  we  do  not  claim  that  the  GP  models  that  we  identify  are 
the  most  accurate  that  can  be  obtained  for  such  a  system,  but  still  serve  as  a  basis  for  comparison 
to  the  EDMD  models  that  are  identified  for  the  same  choice  of  data. 

The  EDMD  procedure  requires  a  selection  of  observables.  Despite  narrowing  down  this  choice 
substantially  by  choosing  to  work  with  linear  and  quadratic  monomials  of  POD  coefficients,  there 
can  still  be  some  additional  ambiguity  that  should  be  explicitly  clarified.  To  begin  with,  one  must 
decide  how  whether  to  subtract  the  mean  flow  from  the  data  before  applying  POD.  This  step  is 
almost  always  performed  when  performing  POD,  partly  because  then  any  reconstruction  of  the 
flow  using  a  linear  combination  of  POD  modes  will  automatically  satisfy  the  required  boundary 
conditions  [66].  Conversely,  it  is  almost  never  done  with  DMD;  doing  so  can  lead  to  an  undesirable 
equivalence  to  taking  a  discrete  Fourier  transform  [27].  Furthermore,  if  one  is  to  subtract  a  “mean”, 
for  this  flow  one  could  conceivable  take  this  to  be  any  of  the  mean  of  the  limit  cycle,  the  mean  of  the 
data,  or  even  the  unstable  equilibrium  velocity  field.  To  emphasize  the  fact  that  we  are  approaching 
this  procedure  from  a  data-driven  perspective,  and  to  be  consistent  with  the  subspace  used  for  both 
procedures,  here  we  first  subtract  the  mean  of  the  data  before  performing  POD  (for  both  GP  and 
EDMD  models).  Note  that  this  is  different  to  what  is  typically  done  when  constructing  GP  models 
for  such  a  system,  where  the  mean  is  most  typically  taken  to  be  the  mean  of  the  limit  cycle  (e.g., 
[95]).  The  mean  and  POD  modes  identified  in  this  manner,  as  well  as  their  relative  energy  content 
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and  the  temporal  evolution  of  their  coefficients,  are  shown  in  Fig.  5.2. 

We  further  note  that  subtracting  the  mean  of  the  velocity  field  before  performing  POD  does  not 
result  in  all  of  the  observables  used  for  EDMD  having  zero  mean,  as  (a,)2  is  always  non-negative. 
When  performing  EDMD,  we  can  also  choose  whether  to  explicitly  include  this  mean  mode  as 
an  observable,  which  amounts  to  allowing  for  a  constant  term  in  Eq.  5.3.6.  We  will  most  often 
choose  to  do  this,  though  will  discuss  this  in  more  detail  later.  On  top  of  all  of  these  details,  the 
regularization  introduced  in  Sect.  5.3.3  introduces  an  additional  parameter  whose  value  must  be 
set,  which  we  will  discuss  in  the  next  section. 

5.4.2  Comparison  between  EDMD  and  GP  models,  regularization,  and  model 
order  dependence 

This  section  analyzes  the  performance  of  both  EDMD  and  GP  models  of  a  range  of  orders  in 
modeling  the  dynamics  of  the  flow  past  a  cylinder.  We  begin  with  some  sample  results.  Fig.  5.3 
shows  the  performance  of  both  Galerkin  projection,  and  the  EDMD  approach  outlined  in  Sect.  5.3, 
in  identifying  a  model  that  can  predict  the  evolution  of  the  first  three  POD  coefficients.  For  this 
system,  the  dynamics  of  these  coefficients  are  known  to  evolve  on  a  paraboloid  [95].  We  observe 
that  the  EDMD  model  is  more  accurate  than  the  GP  model  in  terms  of  obtaining  both  the  correct 
transient  and  limit  cycle  behavior.  Fig.  5.4  shows  that  the  same  findings  hold  when  the  dimension 
of  the  models  increase  to  9.  We  keep  (5  at  a  fixed  value  of  0.5  for  both  of  these  examples,  which 
will  be  the  default  unless  otherwise  mentioned. 

In  order  to  study  more  systematically  the  accuracy  of  models  of  various  order,  we  compare  in 
Fig.  5.5  the  identified  limit  cycle  amplitude  (defined  based  on  the  coefficient  of  the  first  POD  mode) 
and  frequency  for  models  of  order  3-28.  As  well  as  showing  results  for  EDMD  models  identified 
with  /3  =  0.5,  we  additionally  show  the  performance  of  models  identified  with  an  optimized  value  of 
/ 3 ,  where  the  optimal  is  found  based  on  a  direct  search  over  the  range  0  <  /3  <  2.  This  comparison 
shows  that  the  results  are  relatively  insensitive  to  /3,  and  that  the  one  initially  chosen  value  tends 
to  perform  reasonably  well  across  all  model  orders.  It  would  be  possible,  however,  to  develop  more 
sophisticated  methods  to  tune  (3  using  the  data  available.  To  do  this,  one  could  apply  standard 
extrema-seeking  algorithms  to  find  a  value  of  (3  that  allows  the  identified  model  to  best  match  the 
data,  by  some  chosen  metric.  One  could  also  explore  the  possibility  of  replacing  the  / 31  term  in 
Eq.  5.3.10  with  a  diagonal  matrix  with  different  entries,  though  this  then  loses  the  direct  connection 
to  the  joint  least  squares  problem  formulated  in  Eq.  5.3.9. 

In  general,  EDMD  models  are  mode  accurate  than  GP  models  at  predicting  the  limit  cycle 
characteristics  for  models  of  order  less  than  8,  while  GP  models  are  more  accurate  for  high  order, 
except  for  models  of  order  14,  where  the  GP  model  performs  uncharacteristically  poorly.  While 
for  this  simple  system  the  high  order  models  are  not  required  to  obtain  an  accurate  representation 
of  the  system  dynamics,  it  is  important  to  verify  that  the  proposed  algorithm  remains  capable  of 
identifying  stable  and  accurate  models  as  the  model  order  increases,  in  order  to  be  of  use  for  more 
complex  fluids  systems. 

5.4.3  Model  prediction  for  untrained  conditions 

An  important  feature  of  any  identified  model  is  the  ability  to  predict  the  behavior  of  the  system 
along  trajectories  in  phase  space  that  are  not  contained  in  the  data  used  for  model  identification. 
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(k) 


(1) 


Figure  5.2:  (a)  Mean  mode,  (b)  relative  energy  content  and  (c)  cumulative  energy  content  of  the 
first  30  POD  modes  identi^^^|^^^^|^4^ofef(S^lta^.ons  from  near  the  unstable 
equilibrium  to  the  limit  cycle.  (d)-(l)  Contours  of  streamwise  velocity  for  the  first  nine  POD 
modes,  and  their  corresponding  time-varying  coefficients 
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Figure  5.3:  Performance  of  3rd  order  GP  and  EDMD  nonlinear  models  in  predicting  the  evolution 
of  POD  coefficients  for  transitional  flow  past  a  cylinder,  showing  (a)  time  evolution  and  (b)  phase 
portrait  plots.  The  Galerkin  and  EDMD  phase  portrait  models  are  allowed  to  evolve  for  800 
dimensionless  time  units  to  confirm  limit  cycle  behavior 
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Figure  5.4:  Performance  of  9th  order  GP  and  EDMD  nonlinear  models  in  predicting  the  time- 
evolution  of  POD  coefficients  for  transitional  flow  past  a  cylinder,  showing  (a)  time  evolution  and 
(b)  phase  portrait  plots.  Only  the  first  5  modes  are  shown  on  the  left,  while  the  phase  protein 
shows  the  projection  onto  the  first  3  modes.  The  Galerkin  and  eEMD  phase  portrait  models  are 
allowed  to  evolve  for  800  dimensionless  time  units  to  confirm  limit  cycle  behavior 
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Figure  5.5:  Prediction  of  the  limit  cycle  amplitude  (a)  and  frequency  (b)  of  the  first  POD  mode 
from  both  EDMD  and  GP  models 


Figure  5.6:  Performance  of  3rd  order  GP  and  EDMD  nonlinear  models  in  predicting  the  evolution 
of  POD  coefficients  for  transitional  flow  past  a  cylinder,  showing  (a)  time  evolution  and  (b)  phase 
portrait  plots.  The  initial  condition  is  taken  to  be  the  mean  of  the  limit  cycle 


This  ensures  the  practical  utility  of  such  a  model,  and  increases  confidence  that  the  “correct”  system 
dynamics  are  being  captured.  We  show  in  Fig.  5.6  the  performance  of  3rd  order  EDMD  and  GP 
models  in  predicting  the  evolution  of  the  first  3  POD  coefficients,  starting  at  the  mean  of  the  limit 
cycle.  We  observe  that  the  GP  model  more  accurately  captures  the  initial  transient  as  the  flow 
approaches  the  slow  manifold  (paraboloid),  though  as  before  the  EDMD  model  is  more  accurate  in 
predicting  the  subsequent  approach  to  the  limit  cycle.  Note  that  there  are  two  main  reasons  why 
this  model  might  be  somewhat  inaccurate  in  the  untrained  region:  not  only  are  the  dynamics  in 
this  region  untrained,  but  the  POD  basis  is  additionally  no  longer  energetically  optimal  (in  terms 
of  energy)  for  this  dataset,  so  low  order  models  in  particular  might  not  capture  features  that  are 
both  energetically  and  dynamically  important. 
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Figure  5.7:  Performance  of  3rd  order  GP  and  EDMD  nonlinear  models  in  predicting  the  evolution 
of  POD  coefficients  for  transitional  flow  past  a  cylinder,  showing  (a)  time  evolution  and  (b)  phase 
portrait  plots.  The  Galerkin  and  EDMD  phase  portrait  models  are  identified  from  data  that  is 
corrupted  by  noise  of  standard  deviation  a  =  0.05U.  f3  =  0.5 


5.4.4  Noisy  data 

An  important  quality  for  any  modeling  procedure  to  possess  is  a  robustness  to  noisy  data,  such  as 
that  which  might  be  acquired  from  experiments.  In  Fig.  5.7,  we  show  the  results  of  identifying  3rd 
order  models  to  the  same  data  as  considered  previously,  but  corrupted  by  Gaussian  white  noise, 
with  magnitude  of  5%  of  the  freestream  velocity.  It  is  observed  that  the  EDMD  model  retains  it’s 
accuracy,  whereas  the  GP  model  captures  none  of  the  observed  features  of  the  system.  This  result 
is  interesting,  since  the  noise-corrupted  data  is  actually  being  used  to  identify  the  dynamics  in  the 
case  of  EDMD,  whereas  GP  only  uses  the  data  to  obtain  a  spatial  basis.  As  well  as  showing  the 
robustness  of  the  EDMD  modeling  procedure,  these  results  highlight  the  sensitivity  of  GP  to  this 
basis  selection. 

5.4.5  Data  from  a  restricted  spatial  domain 

Another  potential  drawback  of  GP  models  is  that  they  require  spatially  resolved  data  across  a 
wide  domain,  ideally  such  that  the  boundary  conditions  are  constant.  In  Fig.  5.8,  we  attempt  to 
obtain  reduced  order  models  using  a  small  portion  of  the  domain.  Unlike  GP,  EDMD  is  still  able 
to  produce  an  accurate  model  with  such  a  restricted  data  set.  Here,  we  are  only  using  the  data 
collected  from  this  limited  domain  to  compute  POD,  so  the  modes  that  are  computed  will  differ 
somewhat  from  the  modes  shown  in  Fig.  5.2.  The  fact  that  EDMD  models  in  this  case  retain  their 
accuracy  for  both  noisy  data  and  restricted  spatial  domains  suggests  that  this  modeling  framework 
could  be  particularly  useful  for  experimental  data,  which  often  possesses  both  limitations. 

An  advantage  of  EDMD  is  that  it  does  not  require  spatially  resolved  data,  which  is  required 
for  the  computation  of  spatial  derivatives  in  GP.  In  Fig.  5.9,  we  show  that  EDMD  can  produce  an 
accurate  ROM  with  only  data  from  a  small  number  of  sparse,  randomly  chosen  points.  This  gives 
cause  for  optimism  that  the  method  could  additionally  be  used  in  experiments  for  which  only  point 
sensor  measurements  are  available,  rather  than  the  spatially  resolved  flowfields,  which  can  be  more 
costly  and  time  consuming  to  obtain.  While  the  selection  of  the  sensor  locations  in  this  case  was 
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Figure  5.8:  Performance  of  3rd  order  GP  and  EDMD  nonlinear  models  in  predicting  the  evolu¬ 
tion  of  POD  coefficients  for  transitional  flow  past  a  cylinder,  showing  (a)  time  evolution  and  (b) 
phase  portrait  plots.  Only  data  contained  from  inside  the  square  shown  in  (c)  is  used  for  model 
identification.  The  EDMD  models  were  obtained  with  j3  =  0.1 
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Figure  5.9:  Performance  of  3rd  order  GP  and  EDMD  nonlinear  models  in  predicting  the  evolution 
of  POD  coefficients  for  transitional  flow  past  a  cylinder,  showing  (a)time  evolution  and  (b)  phase 
portrait  plots.  Only  data  sampled  at  green  circles  shown  in  (c)  are  used  for  model  identification 


random,  we  observed  empirically  that  best  results  were  obtained  when  the  sensor  locations  were 
spread  across  different  spatial  regions  of  the  domain.  When  using  only  a  small  number  of  point 
sensors,  we  found  that  better  models  were  obtained  when  we  did  not  include  a  constant  term  in 
our  set  of  observables.  We  additionally  did  not  use  any  form  of  regularization. 


5.4.6  Data  from  limited  temporal  sampling 


In  this  section,  we  explore  the  ability  of  EDMD  and  GP  models  in  predicting  the  evolution  of  a 
system  using  only  a  small  amount  of  temporal  data.  As  in  Sect5.4.5,  we  compute  POD  modes 
(including  the  mean  modes)  from  the  subset  of  data  that  is  used.  Note  in  particular  that  this 
means  that  the  POD  basis  that  we  use  as  our  state  space  will  be  different  to  those  identified  from 
the  full  dataset.  In  Fig.  5.10,  we  show  the  results  of  applying  both  EDMD  and  GP  when  using 
only  the  first  40%,  20%,  and  10%  of  the  time  series  of  data  used  for  the  previous  sections.  We  find 
that,  even  when  the  data  available  is  significantly  reduced  and  has  clearly  not  reached  the  limit 
cycle,  the  EDMD  models  still  accurately  predict  the  presence  (and  to  some  extent,  the  location)  of 
the  limit  cycle.  This  is  in  contract  to  the  GP  model,  where  the  limited  data  available  renders  the 
model  qualitatively  incorrect  for  all  amounts  of  data  used.  Note  that  the  results  in  this  section  did 
not  use  any  form  of  regularization. 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


66 


CHAPTER  5.  NONLINEAR  REDUCED-ORDER  MODELS  USING  EDMD 


(e) 


Time 


- Data 

Galerkin  projection  model 
- Extended  DMD  model 


Figure  5.10:  Performance  of  3rd  order  GP  and  eDMD  nonlinear  models  in  predicting  the  evolution 
of  POD  coefficients  for  transitional  flow  past  a  cylinder,  showing  (a)  time  evolution  and  (b)  phase 
portrait  plots.  Models  are  identified  using  only  the  first  400  (top),  200  (center)  and  100  (bottom) 
snapshots  of  data,  as  shown 
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5.5  Discussion  and  conclusions 

This  work  has  introduced  a  method  to  extract  low-order,  nonlinear  models  from  time-resolved  data, 
by  utilizing  an  extension  of  the  DMD  algorithm.  Through  the  use  of  a  variety  of  different  types  of 
data  from  a  simple  example  of  a  nonlinear  fluids  system,  we  have  attempted  to  evaluate  the  relative 
performance  of  this  method  in  a  range,  in  comparison  with  Galerkin  projection. 

We  believe  that  the  advantages  of  the  present  approach  are  numerous.  We  have  demonstrated 
that  this  data-driven  approach  appears  to  be  more  robust  to  noisy  data,  and  exhibits  greater 
flexibility  in  terms  of  the  spatial  extent  and  resolution  of  data  for  which  the  method  produces 
accurate  models.  We  therefore  hold  strong  hopes  that  this  method  will  be  useful  in  particular  to 
model  experimental  data,  which  often  contains  such  limitations. 

Beyond  this,  the  present  method  requires  no  explicit  knowledge  of  the  underlying  governing 
equations,  though  here  knowing  at  least  the  form  of  the  nonlinearity  was  advantageous.  In  par¬ 
ticular,  note  that  the  correct  form  of  the  observables  is  most  likely  important  when  the  nonlinear 
model  is  used  to  extrapolate  a  small  amount  of  data,  and  ultimately  correctly  predict  the  limit 
cycle  in  Sect.  5.4.6.  In  other  words,  if  the  correct  form  of  the  nonlinearity  is  known,  then  even  a 
small  sample  of  data  can  be  sufficient  to  identify  a  model  that  is  at  least  qualitatively  accurate. 
This  finding  has  important  implications  for  more  complicated  systems,  where  sampling  across  all 
possible  regimes  and/or  locations  in  phase  space  may  be  infeasible. 

In  Sect.  5.4.1,  it  was  observed  that  EDMD  models  are  particularly  accurate  when  the  dimension 
of  the  model  is  small,  in  comparison  to  both  GP  models  of  the  same  order,  and  EDMD  models  of 
higher  order.  The  fact  that  GP  models  are  less  accurate  for  lower  model  orders  is  unsurprising,  since 
(without  any  additional  modification)  they  have  no  way  of  accounting  for  the  effect  of  unmodeled 
modes.  EDMD  models,  on  the  other  hand,  can  at  least  account  for  the  time-averaged  effect  of 
modes  that  are  not  explicitly  included  as  variables.  Note,  however,  that  unmodeled  modes  with 
highly  intermittent  dynamics  could  be  impossible  to  accurately  account  for.  As  the  model  order 
increases,  GP  models  should  be  expected  to  become  more  accurate,  as  the  effect  of  unmodeled 
modes  reduces.  Conversely,  we  find  that  EDMD  models  can  become  less  accurate  and  robust, 
possibly  due  to  the  rapid  increase  in  the  number  of  parameters  requiring  to  be  fit.  Note  that  a 
phenomena  has  been  observed  previously  in  polynomial  identification  [102]. 

The  mathematical  machinery  behind  DMD  and  it’s  variants  involves  a  least  squares  (or  mini¬ 
mum  norm,  in  the  underconstrained  case)  fit  of  data.  In  this  sense,  the  method  presented  in  this 
work  shares  strong  similarities  with  a  number  of  other  methods  that  have  been  used  to  identify 
nonlinear  systems  in  fluids,  such  as  polynomial  identification  [102],  modified  quadratic  stochastic  es¬ 
timation  [94]  Volterra  system  identification  [11],  and  linear  parameter  varying  system  identification 
[60],  which  all  employ  a  similar  type  of  least-squares  coefficient-fitting  methodology  in  their  respec¬ 
tive  algorithms.  The  recently  developed  SINDy  [24]  (sparse  identification  of  nonlinear  dynamics) 
also  employs  a  similar  basic  framework,  though  with  an  additional  sparsity-promoting  algorithm 
to  decide  upon  the  appropriate  observables  from  a  larger  library,  which  is  particularly  useful  when 
the  form  of  the  nonlinearity  is  unknown.  In  this  sense,  we  do  not  claim  that  the  algorithm  we  use 
is  exceptionally  novel,  but  hope  that  this  exposition,  both  in  terms  of  the  connections  of  DMD 
and  Koopman  theory,  and  in  highlighting  a  number  of  situations  where  such  a  modeling  approach 
might  be  advantageous  or  necessary.  On  the  former  point,  one  could  hope  that  EDMD,  with  the 
correct  choice  of  observables,  can  allow  for  accurate  prediction  not  just  of  Koopman  eigenvalues 
and  eigenmodes,  but  also  Koopman  eigenfunctions.  There  is  indeed  strong  evidence  to  suggest  that 
identified  Koopman  eigenfunctions  can  indeed  give  an  accurate  characterization  of  the  dynamics  of 
the  cylinder  wake  near  the  limit  cycle  [7].  We  do  not  make  any  claims  that  we  have  achieved  such 
convergence  here  over  the  full  transient  regime,  but  rather  show  that  the  results  can  still  be  useful 
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from  a  practical,  model-reduction  standpoint.  We  do  not  make  any  claims  that  we  have  achieved 
such  convergence  here,  but  show  that  the  results  can  still  be  useful  from  a  practical,  model-reduction 
standpoint.  As  well  as  pursuing  this  direction,  further  work  could  seek  to  utilize  such  models  for 
purposes  of  flow  control.  Preliminary  investigations  also  suggest  that  such  EDMD-based  models 
can  be  effective  at  identifying  stability  properties,  such  as  eigenvalues  of  the  linearized  system  at 
equilibrium  points. 
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Chapter  6 

Characterizing  and  correcting  for  the 
effect  of  sensor  noise  in  the  dynamic 
mode  decomposition 


Dynamic  mode  decomposition  (DMD)  provides  a  practical  means  of  extracting  insightful  dynamical 
information  from  fluids  datasets.  Like  any  data  processing  technique,  DMD’s  usefulness  is  limited 
by  its  ability  to  extract  real  and  accurate  dynamical  features  from  noise-corrupted  data.  Here  we 
show  analytically  that  DMD  is  biased  to  sensor  noise,  and  quantify  how  this  bias  depends  on  the  size 
and  noise  level  of  the  data.  We  present  three  modifications  to  DMD  that  can  be  used  to  remove  this 
bias:  (i)  a  direct  correction  of  the  identified  bias  using  known  noise  properties,  (ii)  combining  the 
results  of  performing  DMD  forwards  and  backwards  in  time,  and  (iii)  a  total  least-squares-inspired 
algorithm.  We  discuss  the  relative  merits  of  each  algorithm,  and  demonstrate  the  performance 
of  these  modifications  on  a  range  of  synthetic,  numerical,  and  experimental  datasets.  We  further 
compare  our  modified  DMD  algorithms  with  other  variants  proposed  in  recent  literature. 

6.1  Introduction 

With  advances  in  both  experimental  techniques  and  equipment,  and  computational  power  and 
storage  capacity,  researchers  in  fluid  dynamics  can  now  generate  more  high-fidelity  data  than  ever 
before.  The  presence  of  increasingly  large  data  sets  calls  for  appropriate  data  analysis  techniques, 
that  are  able  to  extract  tractable  and  physically  relevant  information  from  the  data.  Dynamic  mode 
decomposition  allows  for  the  identification  and  analysis  of  dynamical  features  of  time-evolving  fluid 
flows,  using  data  obtained  from  either  experiments  or  simulations.  In  contrast  to  other  data- 
driven  modal  decompositions  such  as  the  proper  orthogonal  decomposition  (POD),  DMD  allows 
for  spatial  modes  to  be  identified  that  can  be  directly  associated  with  characteristic  frequencies  and 
growth/decay  rates.  Following  its  conception,  DMD  was  quickly  shown  to  be  useful  in  extracting 
dynamical  features  in  both  experimental  and  numerical  data  [125,  126].  It  has  subsequently  been 
used  to  gain  dynamic  insight  on  a  wide  range  of  problems  arising  in  fluid  mechanics  [e.g.,  128,  127, 
70]  and  other  fields  [e.g.,  59]. 

DMD  has  a  strong  connection  to  Koopman  operator  theory  [81,  88],  as  exposed  in  [119],  and 
further  reviewed  in  [89],  which  can  justify  its  use  in  analyzing  nonlinear  dynamical  systems.  Since 
its  original  formulation,  numerous  modifications  and  extensions  have  been  made  to  DMD.  [27] 
highlights  the  connection  that  DMD  shares  with  traditional  Fourier  analysis,  as  well  as  proposing 
an  optimized  algorithm  that  recasts  DMD  as  an  optimal  dimensionality  reduction  problem.  This 


69 

DISTRIBUTION  A:  Distribution  approved  for  public  release. 


70 


CHAPTER  6.  CORRECTING  FOR  SENSOR  NOISE  IN  DMD 


concept  of  finding  only  the  dynamically  important  modes  has  also  been  considered  in  subsequent 
works  of  [164]  and  [71].  All  of  these  works  are  motivated,  in  part,  by  the  fact  that  by  default  DMD 
will  output  as  many  modes  as  there  are  pairs  of  snapshots  (assuming  that  the  length  of  the  snapshot 
vector  is  greater  than  the  number  of  snapshots),  which  is  arbitrary  with  respect  to  the  dynamical 
system  under  consideration.  In  reality,  one  would  prefer  to  output  only  the  modes  and  eigenvalues 
that  are  present  (or  dominant)  in  the  data.  When  the  data  is  corrupted  by  noise  (as  will  always  be 
the  case  to  some  degree,  especially  for  experimental  data),  this  process  becomes  nontrivial,  since 
noisy  data  might  have  a  numerical  rank  far  larger  than  the  dimension  of  the  governing  dynamics  of 
the  system.  Further  to  this,  one  cannot  expect  to  have  a  clean  partition  into  modes  that  identify 
true  dynamical  features,  and  those  which  consist  largely  of  noise. 

Simple  ways  of  achieving  this  objective  can  involve  either  first  projecting  the  data  onto  a  smaller 
dimensional  basis  (such  as  the  most  energetic  POD  modes)  before  applying  DMD,  or  by  choosing 
only  the  most  dynamically  important  DMD  modes  after  applying  DMD  to  the  full  data.  One 
can  also  truncate  the  data  to  a  dimension  larger  than  the  assumed  dimension  of  the  dynamics, 
and  then  apply  a  balanced  truncation  to  the  resulting  dynamical  system  to  obtain  the  desired 
reduced  order  model.  This  approach  is  sometimes  referred  to  as  over-specification  in  the  system 
identification  community  [see,  e.g.,  75].  Keeping  a  higher  dimension  of  data  than  that  of  the  assumed 
dynamics  can  be  particularly  important  for  input-output  systems  that  have  highly  energetic  modes 
that  are  not  strongly  observable  or  controllable  [112].  Ideally,  any  algorithm  that  restricts  the 
number  of  DMD  modes  that  are  computed  should  also  additionally  be  computationally  efficient. 
A  fast  method  to  perform  DMD  in  real  time  on  large  datasets  was  recently  proposed  in  [62],  while 
a  library  for  efficient  parallel  implementation  of  number  of  common  modal  decomposition  and 
system  identification  techniques  is  described  in  [16].  An  extension  of  DMD  that  potentially  allows 
for  better  representation  of  nonlinear  data  has  also  recently  been  proposed  [160],  and  although 
the  computational  costs  increase  dramatically  with  the  dimension  of  the  system,  a  kernel  method 
described  in  [161]  reduces  the  cost  to  be  comparable  to  standard  DMD. 

One  of  the  major  advantages  of  DMD  over  techniques  such  as  global  stability  analysis,  for  exam¬ 
ple,  is  that  it  can  be  applied  directly  to  data,  without  the  need  for  the  knowledge  or  construction  of 
the  system  matrix,  which  is  typically  not  available  for  experiments  [126].  For  this  reason,  analysis 
of  the  sensitivity  of  DMD  to  the  type  of  noise  typically  found  in  experimental  results  is  of  particular 
importance.  The  effects  of  noise  on  the  accuracy  of  the  DMD  procedure  was  systematically  inves¬ 
tigated  in  the  empirical  study  of  [44],  for  the  case  of  a  synthetic  waveform  inspired  by  canonical 
periodic  shear  flow  instabilities.  More  recently,  [101]  have  extended  this  type  of  analysis  to  more 
complex  data  with  multiple  frequencies,  as  might  be  found  in  typical  fluids  systems.  The  present 
work  builds  upon  these  previous  studies  by  analytically  deriving  an  expression  that  explicitly  shows 
how  DMD  should  be  affected  by  noise,  for  the  case  where  the  noise  is  assumed  to  be  sensor  noise 
that  is  uncorrelated  with  the  dynamics  of  the  system.  Our  analysis  complements  the  “noise-robust” 
DMD  formulation  in  [63]  by  explicitly  quantifying  the  influence  of  noise  on  DMD.  Further,  while 
our  analysis  is  consistent  with  the  total  least-squares  formulation  in  [63],  we  use  the  insights  gained 
from  our  analysis  to  develop  alternative  techniques  to  total  least-squares  DMD  that  may  be  prefer¬ 
able  in  certain  applications.  Ultimately,  the  availability  of  multiple  “noise-aware”  DMD  algorithms 
allows  the  user  to  approach  dynamical  analysis  of  noisy  data  from  multiple  angles,  thus  garnering 
more  confidence  in  the  computations.  We  note  that  the  case  of  process  noise,  where  noise  can 
interact  with  the  dynamics  of  the  system,  is  also  the  subject  of  recent  work  [8]. 

Our  analysis  uses  a  recent  characterization  of  DMD  [149],  which  highlights  the  connection  of 
DMD  to  related  techniques  that  are  used  in  other  communities  for  the  extraction  of  dynamical 
information  from  data.  Many  linear  system  identification  techniques  are  closely  related  in  that 
they  are  based  around  singular  value  decomposition  of  a  data  matrix;  aside  from  DMD  there  is 
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the  eigensystem  realization  algorithm  [72]  and  balanced  proper  orthogonal  decomposition  [112],  for 
example.  Indeed,  the  origin  of  such  an  approach  seems  to  date  back  to  the  work  of  [65]. 

In  this  work,  we  first  show  that  the  dominant  effect  of  noise  on  DMD  is  often  deterministic. 
This  not  only  allows  us  to  accurately  predict  its  effect,  but  also  allows  for  a  correction  to  be 
implemented  to  recover  the  noise-free  dynamics.  As  well  as  directly  correcting  for  the  noise,  we 
present  two  other  modifications  of  DMD,  that  both  are  able  to  remove  this  bias  without  needing  to 
know  the  noise  characteristics.  Sect.  6.2  develops  the  theory  that  characterizes  the  effect  of  noise 
on  DMD,  which  subsequently  motivates  the  formulation  of  our  modified  algorithms,  which  we  term 
noise-corrected  DMD  (ncDMD),  forward-backward  DMD  (fbDMD)  and  total  least-squares  DMD 
(tlsDMD).  In  Sect.  6.3,  we  analyze  the  performance  of  these  algorithms  on  a  number  of  synthetic 
data  sets,  which  are  corrupted  by  Gaussian  white  noise.  We  additionally  investigate  how  the 
algorithms  perform  on  data  with  both  sensor  and  process  noise.  In  Sect.  6.4,  we  use  numerical  and 
experimental  data  from  flow  past  a  cylinder  undergoing  periodic  vortex  shedding,  to  demonstrate 
the  utility  of  the  proposed  modifications  of  DMD  for  real  fluids  data. 

6.2  Characterizing  noise  in  dynamic  mode  decomposition 

This  section  details  the  methodology  that  is  used  to  analyze  the  effect  of  noise  in  DMD.  After 
introducing  DMD  in  Sect.  6.2.1,  the  effect  of  sensor  noise  in  the  data  on  the  results  of  DMD  is 
studied  in  Sect.  6.2.2,  which  in  particular  shows  that  DMD  is  biased  to  sensor  noise.  Sects.  6.2.3- 
6.2.5  formulate  three  different  modifications  of  the  DMD  algorithm  that  are  designed  to  remove 
this  bias. 

6.2.1  Dynamic  mode  decomposition 

DMD  has  undergone  a  number  of  formulations,  interpretations  and  modifications  since  its  inception. 
Common  to  all  methods  is  the  requisite  collection  and  arrangement  of  data,  summarized  now. 
Suppose  we  collect  snapshots  of  data  x*,  which  we  assemble  as  columns  in  the  data  matrix  Z. 
For  fluids  systems  x*  will  typically  be  a  velocity  field  snapshot,  but  more  generally  it  is  a  vector 
of  observations  of  an  evolving  dynamical  system  at  a  given  time.  From  Z ,  we  select  all  pairs  of 
columns  that  are  sampled  at  a  time  difference  At  apart,  and  place  them  into  the  matrices  X  and  Y 
(where  the  data  in  a  given  column  of  Y  was  collected  At  after  the  equivalent  column  of  A).  Note 
that  if  Z  consists  of  a  sequential  time-series  of  data,  then  X  and  Y  are  simply  Z  with  the  last  and 
first  columns  excluded,  respectively.  Let  X  and  Y  each  be  n  by  m  matrices,  so  we  have  m  pairs  of 
snapshots,  each  of  size  n.  By  not  explicitly  requiring  a  single  time-series  of  data,  we  allow  for  larger 
or  irregular  time  gaps  between  snapshot  pairs,  the  concatenation  of  data  from  multiple  time-series, 
and  for  the  removal  of  any  corrupted  or  spurious  data.  Recently,  [149]  proposed  an  interpretation 
of  DMD  modes  and  eigenvalues  as  the  eigendecomposition  of  the  matrix 


A  =  YX+,  (6.2.1) 

where  X+  denotes  the  Moore-Penrose  pseudoinverse  of  a  matrix  X.  While  this  is  a  succinct 
interpretation,  and  one  which  will  be  useful  in  the  ensuing  analysis,  it  is  typically  not  an  efficient 
(or  even  feasible)  means  of  performing  DMD  (as  discussed  in  [149]).  This  is  especially  true  when 
n  3>  m,  which  is  often  the  case  in  high-dimensional  fluids  systems.  Instead,  since  X  and  Y  have 
rank  at  most  min(m,  n),  it  is  typically  more  efficient  to  first  project  the  data  onto  a  subspace  that 
is  (at  most)  of  this  dimension.  One  way  to  do  this  is  by  projecting  the  original  snapshots  onto 
the  POD  modes  of  the  data,  which  is  implicitly  done  in  all  DMD  algorithms.  Note  that  the  POD 
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modes  of  X  are  the  columns  of  U  in  the  singular  value  decomposition  X  =  UYV*  (though  typically 
POD  is  performed  after  first  subtracting  the  temporal  mean  of  the  data,  which  is  not  done  here). 
We  present  here  a  typical  algorithm  to  compute  DMD,  that  is  most  similar  to  that  proposed  in 
[149]  as  exact  DMD. 

Algorithm  2  (DMD). 

1.  Take  the  reduced  singular  value  decomposition  (SVD)  of  X ,  letting  X  =  UT,V* . 

2.  (Optional)  Truncate  the  SVD  by  only  considering  the  first  r  columns  of  U  and  V,  and  the 
first  r  rows  and  columns  of  £  (with  the  singular  values  ordered  by  size),  to  obtain  Ur,  £r  , 
and  Vr 

3.  Let  A  :=  UfYVrE-1 

4 ■  Find  the  eigenvalues  )Jn  and  eigenvectors  Wi  of  A,  with  Awt  =  HiWi, 

5.  Every  nonzero  m  is  a  DMD  eigenvector,  with  a  corresponding  DMD  mode  given  by  ipi  := 
H^YVrE^Wi. 

This  method  is  similar  to  the  original  formulation  in  [126],  but  for  the  fact  that  in  step  5 
the  DMD  modes  are  no  longer  restricted  to  lie  within  the  column  space  of  X.  We  also  explicitly 
provide  the  optional  step  of  truncating  the  SVD  of  X,  which  might  be  done  if  the  system  is  known 
to  exhibit  low  dimensional  dynamics,  or  in  an  attempt  to  eliminate  POD  modes  that  contain 
only  noise.  We  note  that  this  is  not  the  only  means  to  reduce  the  dimension  of  the  identified 
system  dynamics,  nor  is  it  necessarily  optimal.  Indeed,  [164]  develops  a  method  that  optimizes  the 
projection  basis  in  parallel  while  performing  a  DMD-like  eigendecomposition.  [71]  takes  a  different 
approach,  seeking  a  small  number  of  nonzero  modes  from  the  full  eigendecomposition  that  best 
approximate  the  system  dynamics.  An  empirical  comparison  between  these  various  dimensionality- 
reduction  techniques  will  be  given  in  Sect  6.3.3.  Note  that  the  continuous  eigenvalues  A cj  of  the 
system  are  related  to  the  discrete  time  eigenvalues  identified  via  DMD  via  Acj  =  log(ni)/ At.  The 
growth  rate  7 j  and  frequency  07  associated  with  DMD  mode  ipi  are  then  given  by  A d  =  7 ,  +  iuoi. 

The  matrix  A  is  related  to  A  in  Eq.  (6.2.1)  by  A  =  U(AUr.  While  A  can  be  viewed  as  an 
approximating  linear  propagation  matrix  in  Rn  (i.e.,  the  space  of  original  data  vectors),  A  is  the 
equivalent  propagation  matrix  in  the  space  of  POD  coefficients,  which  we  will  sometimes  refer  to 
as  POD  space.  Another  interpretation  of  A  is  that  it  is  the  spatial  correlation  matrix  between  the 
POD  modes  Ur,  and  the  same  POD  modes  shifted  by  the  assumed  dynamics  A  [126].  If  we  let 
Xfc  =  be  the  representation  of  a  given  snapshot  x  in  the  POD  basis  and  let  X  =  UfX  and 

Y  =  U*Y,  then  it  is  easy  to  verify  that  the  equivalent  of  Eq.  (6.2.1)  in  POD  space  is 


A  =  YX+. 


(6.2.2) 


Eq.  (6.2.2)  will  be  useful  for  the  subsequent  analysis  performed  in  this  paper. 

6.2.2  Sensor  noise  in  DMD 

In  this  work  we  use  the  term  sensor  noise  to  describe  additive  noise  that  affects  only  our  measure¬ 
ments  of  a  given  system,  and  does  not  interact  with  the  true  dynamics.  If  we  have  a  discrete-time 
dynamical  system 

x(f  +  At)  =  F[x(t)], 
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then  we  assume  that  our  measurements  take  the  form 

xm(t)  =  x(t)  +  n  (t), 

where  n(t)  is  a  random  noise  vector.  For  the  purposes  of  this  paper,  we  will  take  each  component  of 
n(f)  to  be  independent  and  normally  distributed  with  zero  mean  and  a  given  variance.  For  instance, 
zero  mean.  I  don’t  think  we  need  to  be  pedantic  about  the  probability  theory  here,  but  probably 
good  to  be  a  little  more  specific  than  this.  With  X  and  Y  as  described  in  Sect.  6.2.1,  suppose 
that  we  measure  Xm  =  X  +  Nx  and  Ym  =  Y  +  Ny,  where  Nx  and  Ny  are  random  matrices  of 
sensor  noise.  Note  that  some  (or  most)  columns  of  random  data  in  Nx  might  also  be  in  Ny,  but 
shifted  to  a  different  column.  We  will  assume  that  the  noise  is  independent  of  the  true  data,  and 
is  independent  in  both  space  and  time,  so  that  each  element  of  a  given  noise  matrix  is  a  random 
variable  taken  from  a  fixed  zero- mean  normal  distribution.  From  Eq.  (6.2.2),  the  measured  DMD 
matrix  Am  can  be  computed  from 

Am  =  Ymx+  =  (Y  +  Ny)(x  +  Nx)+ 

=  (Y  +  Ny)(x  +  Nx)*[(X  +  Nx)(X  +  Nx)*}+ 

=  (YX*  +  NyX*  +  YN*X  +  NyN*x)  [ll*  +  NxX*  +  XN*x  +  NXNX 1 +  ,  (6.2.3) 

where  we  have  used  the  identity  M+  =  M*(MM*)+ .  Note  that  here  the  r  notation  means  that  the 
data  is  expressed  in  the  POD  basis  obtained  from  the  noisy  data.  We  perform  our  analysis  in  this 
POD  space  rather  than  with  the  original  data  to  allow  for  truncation  of  low  energy  modes,  and 
because  the  computation  of  the  pseudoinverse  X+  can  be  prohibitive  for  large  datasets.  We  expect 
that  the  presence  of  noise  should  result  in  some  error  in  the  computation  of  Am  (in  comparison  to 
the  noise  free  matrix  A)  and  thus  some  amount  of  error  in  the  computed  DMD  eigenvalues  and 
modes.  Since  elements  of  Am  are  statistical  quantities  dependent  on  the  noise,  it  will  make  sense  to 
compute  statistical  properties  of  the  matrix.  We  begin  by  computing  E[Am],  the  expected  value  of 
the  computed  DMD  matrix.  Provided  that  we  have  truncated  any  POD  modes  with  zero  energy, 
XX*  should  be  invertible.  If  the  noise  terms  are  sufficiently  small,  then  we  can  make  use  of  the 
matrix  perturbed  inverse  expansion  (M  +  P)^1  =  AM1  —  M~lPM~ 1  +  . . . ,  where  higher  order 
terms  will  be  small  for  M  3>  P.  Eq.  6.2.3  then  becomes 

Am  =  (YX*  +  NyX*  +  YN*X  +  NyN*x)(XX*)-1  \l  -  (NxX*  +  XN*x  +  NxN*x)(XX*)~l  +  ...  . 

(6.2.4) 

Taking  the  expected  value  of  Eq.  (6.2.4),  we  may  classify  the  terms  into  three  categories:  a 
deterministic  terms  that  does  not  involve  Nx  or  Ny  (which  ends  up  being  A),  terms  involving  one 
or  three  noise  matrices,  which  will  have  expected  values  of  0  (e.g.,  NyX* (XX*)~l),  and  terms  which 
involve  two  or  four  noise  matrices.  It  is  terms  this  latter  category  that  may  have  non-zero  expected 
values,  and  thus  bias  the  result  of  applying  DMD  to  noisy  data.  Discarding  terms  containing  a 
single  noise  matrix,  and  additionally  discarding  higher  order  terms  from  the  expansion,  we  have 

E (Am)  =  A(I  -  EiNxN^iXX*)-1)  +  E (NyX+Nx)X+  +  E (NyX+XN^XX*)-1 
+  YE(N*x(XX*)~lNx)X+  +  YE(N*x(XX*)~lXN*x(XX*)~l) 

+  E  \NyNx(XX*)—^(I  -  NxNxiXX*)-1)  ,  (6.2.5) 

where  we  have  noted  that  YX* (XX*)^1  =  YX+  =  A.  Assuming  that  the  noise  is  sufficiently 
small  compared  with  the  true  data,  we  can  further  neglect  the  term  involving  four  noise  matrices. 
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Figure  6.1:  Illustrative  diagram  showing  how  the  error  in  estimation  of  a  given  quantity  can  be 
decomposed  into  bias  error  (being  the  difference  between  the  true  and  expected  value  of  the  iden¬ 
tified  quantity),  and  random  error  (representing  the  fluctuation  in  the  estimated  quantity  between 
different  noise  realizations) 


The  largest  of  the  remaining  terms  will  be  that  which  contains  the  product  NxN^-.  The  remaining 
terms  do  not  necessarily  have  zero  mean,  but  for  the  purposes  of  this  investigation  will  be  neglected. 
Our  results  will  demonstrate  that  this  simplification  is  justifiable.  This  reduces  Eq.  (6.2.5)  to  the 
following  expression,  relating  the  identified  and  true  DMD  matrices: 

E(Am)  =  A(I  -  E (NxN^XX*)-1).  (6.2.6) 

It  might  seem  surprising  that  Eq.  (6.2.6)  contains  Nx,  but  not  Ny ■  The  reason  for  this  will 
become  apparent  in  Sect.  6.2.5,  where  casting  DMD  in  an  optimization  framework  shows  that  the 
standard  algorithm  is  optimal  only  when  assuming  that  all  of  the  noise  is  in  Y,  but  not  X.  From 
a  mathematical  point  of  view,  it  is  because  the  expression  A  =  YX+  is  linear  in  Y  but  not  in  X, 
which  is  why  perturbations  to  X  do  not  have  to  propagate  through  the  equation  in  an  unbiased 
manner.  Note  that  the  same  analysis  can  be  performed  without  transforming  into  POD  space  (i.e., 
without  the  7  notation),  with  the  analogous  expression  to  Eq.  (6.2.6)  being 

E (An)  =  A(I  -  E (NxN^iXX*)-1),  (6.2.7) 

subject  to  XX*  being  invertible.  For  systems  where  the  size  of  each  snapshot  is  larger  than  the 
number  of  snapshots  (i.e,  n  >  m,  which  is  typical  for  fluids  systems),  XX*  will  not  be  invertible, 
thus  motivating  our  choice  to  work  in  POD  space.  Moreover,  one  might  want  the  option  to  truncate 
all  but  a  certain  number  of  POD  modes,  in  order  to  obtain  a  low-dimensional  model  for  the  dominant 
system  dynamics.  Up  until  this  point,  we  have  not  made  a  distinction  between  the  POD  modes  of 
the  clean  data,  X ,  and  the  noisy  measured  data,  Xm.  with  the  latter  typically  being  all  that  we 
have  access  to.  This  issue  will  be  explicitly  addressed  in  Sect.  6.2.3. 

Eq.  (6.2.6)  shows  that  DMD  is  biased  to  sensor  noise.  In  practice,  the  importance  of  this 
finding  will  depend  on  how  the  magnitude  of  this  bias  compares  to  the  random  component  of 
error,  that  will  fluctuate  with  different  samples  of  random  noise.  Figure  6.1  shows  an  illustration 
of  how  bias  and  random  components  of  error  contribute  to  the  total  error  in  the  estimation  of 
some  quantity  from  noisy  data.  Appendix  1  provides  scaling  arguments  that  suggest  that  the  bias 
will  be  the  dominant  component  of  error  in  DMD  whenever  m1/2 SN R  >  n1/2,  where  SNR  is  the 
signal-to-noise  ratio.  When  this  is  the  case,  it  would  be  particularly  advantageous  if  one  had  access 
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to  a  bias-free  alternative  to  DMD.  The  remainder  of  this  section  will  present  a  number  of  such 
alternatives. 

6.2.3  Direct  correction  of  sensor  noise  bias  in  DMD 

Referring  back  to  Eq.  (6.2.6),  we  can  form  a  bias-free  estimate  of  the  true  DMD  matrix  A  via 

A  «  Am(I  -  E (NxNZXXX*)-1)-1.  (6.2.8) 

Making  this  modification  in  practice  requires  an  accurate  estimate  of  both  the  noise  covariance, 
E(NxN^),  and  the  true  data  covariance,  XX*,  in  POD  space.  For  noise  that  is  sufficiently  small, 
we  can  utilize  the  approximation 

XX*  =  U*XX*U  «  U*mXmX*mUm  =  E2m,  (6.2.9) 

where  UmEmV^  is  the  singular  value  decomposition  of  the  noisy  data,  Xm.  This  allows  for  us 
to  express  the  bias  of  DMD  in  terms  of  quantities  that  are  measurable  from  noisy  data.  The 
assumption  that  XX*  =  (Xm  —  Nx)(Xm  —  Nx)*  ~  XrnX*n  can  be  further  refined  by  retaining  the 
NxN^  term,  but  for  small  noise  this  higher  order  term  will  typically  be  small  enough  to  neglect 
after  being  inserted  into  Eq.  (6.2.8).  The  assumption  that  U  ~  Urn  will  largely  be  justified  by 
means  of  results  that  show  the  utility  of  this  analysis.  Analyzing  the  precise  relationship  between 
U  and  Um  in  more  detail  is  beyond  the  scope  of  this  work,  and  is  indeed  an  active  area  of  research. 
We  direct  the  interested  reader  to  relevant  results  in  perturbed  SVD’s  [137,  138,  9],  random  inner 
product  matrices  [28,  141],  and  POD-type  operations  on  noisy  data  [48,  133,  169]. 

Assuming  that  the  noise  is  uniform  as  well  as  spatially  and  temporally  independent,  then 
E(NxNx)  =  E(U*NxN^-U)  =  U*mo2NU  =  ma2NI,  where  a2N  is  the  variance  of  each  independent 
component  of  the  noise  matrix.  With  this  assumption,  and  the  approximation  given  in  Eq.  (6.2.9), 
Eq.  (6.2.8)  becomes 

A  ~  Am(I  —  mcr^E”2)-1.  (6.2.10) 

If  the  noise  is  sufficiently  small,  then  a  perturbed  inverse  approximation  gives 

A  ~  Am(I  +  mfj^rEm2).  (6.2.11) 

We  thus  have  derived  a  correction  to  the  bias  that  is  present  in  the  original  DMD  matrix  Am  due 
to  the  effect  of  sensor  noise.  We  note  that  this  approximation  relies  on  an  accurate  knowledge  of 
the  noise  covariance  matrix.  There  are  numerous  means  to  estimate  noise  properties  from  data,  see 
[108],  for  example.  The  approximations  used  in  deriving  this  expression  also  rely  on  the  magnitude 
of  the  noise  being  smaller  than  that  of  the  true  data  within  each  non-truncated  POD  mode.  We 
now  state  explicitly  the  algorithm  by  which  we  can  correct  for  the  effect  of  sensor  noise  in  the  DMD 
algorithm,  which  we  refer  to  as  noise-corrected  DMD,  or  ncDMD: 

Algorithm  3  (Noise-corrected  DMD  (ncDMD)). 

1.  Compute  Am  from  the  measured  data  as  per  steps  1-3  of  Algorithm  2 

2.  Compute  the  approximation  of  A  from  Eq.  (6.2.10) 

3.  Compute  the  DMD  eigenvalues  and  modes  via  steps  4-5  of  Algorithm  2,  using  the  bias-free 
estimate  of  A. 
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As  was  also  noted  in  Sect  6.2.2,  we  could  have  performed  all  of  the  above  analysis  without 
first  projecting  onto  the  space  of  POD  coefficients,  which  gives  us  the  following  as  analogous  to 
Eqs.  (6.2.10)  and  (6.2.11)  respectively,  subject  to  the  appropriate  inverses  existing: 

A  ps  Am(I  -  maUXmXfJ-1)-1  «  A,m(I  +  ma%(XmXfn)~1)  ■  (6.2.12) 

While  this  approach  might  be  computationally  prohibitive  for  many  applications  of  DMD  (since  it 
requires  inverting  large  n  x  n  matrices),  it  could  in  theory  be  more  accurate,  since  it  doesn’t  rely 
on  any  assumption  that  the  POD  modes  for  the  measured  and  true  data  are  sufficiently  close  to 
each  other.  Note  again  that  XX*  can  only  be  invertible  if  m  >  n,  as  otherwise  it  cannot  be  full 
rank. 

6.2.4  Forward- backward  DMD 

If  we  were  to  swap  the  data  in  X  and  Y.  then  (for  suitably  well  behaved  data)  we  should  expect  to 
identify  the  inverse  dynamical  system,  with  state  propagation  matrix  Bm  (or  Bm  in  POD  space), 
which  approximates  the  true  dynamics  B  (and  B).  Note  that  it  is  not  guaranteed  that  the  dynamics 
of  the  original  system  are  invertible,  but  this  assumption  should  not  be  too  restrictive  for  the 
majority  of  physical  systems  under  consideration  (particularly  after  projection  onto  an  appropriate 
POD  subspace).  It  is  argued  in  Appendix  1  that  sensor  noise  has  the  effect  of  shifting  the  computed 
DMD  eigenvalues  to  appear  to  be  more  stable  than  they  actually  are  (i.e.,  moving  them  further 
inside  the  unit  circle).  Since  our  analysis  was  independent  of  the  nature  of  the  data,  we  should 
expect  the  same  effect  to  be  present  for  the  computation  of  the  inverse  system.  However,  if  B  is 
invertible,  then  we  should  have  B  =  A-1,  meaning  that  we  should  be  able  to  compute  an  estimate  of 
the  forward-time  propagation  matrix  using  backward-time  DMD,  via  A^ck  =  I?”1 .  However,  given 
that  the  eigenvalues  of  Bm  should  have  their  growth  rates  underestimated,  those  of  the  eigenvalues 
of  A^ck  will  be  overestimated.  Specifically,  from  consideration  of  Eq.  (6.2.6),  we  have 

E (Bm)  *b(i-  E (NxN*x)(XX*)-^  , 

and  so 

AbZk  ~  (i  ~  E {NxN*x){XX*)-^j  1  A,  (6.2.13) 

where  we  are  using  the  fact  that  the  noise  and  POD  energy  components  are  the  same  for  forward- 
and  backward-DMD.  We  can  then  combine  estimates  of  the  dynamics  from  forward-  and  backward¬ 
time  DMD  to  obtain 

AmAbZk  =  A(l-  E  {NxN*x){XX*)-^j  (i  -  E  (NxN*x){XX*)-1^  1  A  =  A2.  (6.2.14) 

We  thus  have  the  estimate 

A^{AmAbZk)1/2.  (6.2.15) 

Note  that  this  square  root  will  in  general  be  non-unique,  and  thus  determining  which  root  is  the 
relevant  solution  could  be  nontrivial.  One  reasonable  method,  if  there  is  any  ambiguity,  is  to 
take  the  square  root  which  is  closest  to  Am  (or  AbZk )•  See  [53]  for  a  more  detailed  discussion 
of  the  computation  of  matrix  square  roots.  As  an  aside,  note  that  if  we  assume  that  we  know 
the  equivalent  continuous  time  matrices  Afm  =  log(Am)/At  and  =  log(A^cfc)/Af,  then  the 

equivalent  of  Eq.  (6.2.15)  is 

I  2c,back\ 

^  ~  2  '  -™-m  )• 

We  are  now  in  a  position  to  formalize  this  algorithm,  which  we  refer  to  as  forward-backward,  DMD 
or  fbDMD. 
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Algorithm  4  (forward-backward  DMD  (fbDMD)). 

1.  Compute  Am  from  the  measured  data  as  per  steps  1-3  of  Algorithm  2 

2.  Compute  Bm  from  the  measured  data  as  per  steps  1-3  of  Algorithm  2,  where  X  and  Y  are 
interchanged 

3.  Compute  the  approximation  of  A  from  Eq.  (6.2.15) 

4-  Compute  the  DMD  eigenvalues  and  modes  via  steps  4-5  of  Algorithm  2,  using  the  improved 
estimate  of  A  from  step  f. 

Note  that  in  the  case  where  most  data  snapshots  are  in  both  X  and  Y  (e.g.,  for  a  sequential 
time  series  of  data)  we  can  reduce  the  computational  cost  of  steps  1-2  in  Algorithm  4  by  first  taking 
the  SVD  of  the  entire  data  set,  and  then  working  in  the  space  of  the  resulting  POD  modes. 

6.2.5  Total  least-squares  DMD 

For  the  case  where  the  number  of  snapshots,  m,  is  greater  than  the  size  of  each  snapshot,  n, 
the  DMD  matrix  A  can  be  interpreted  as  the  least-squares  solution  to  the  overdetermined  system 
AX  =  Y.  When  n  >  m,  then  the  solution  for  the  now  underdetermined  system  is  the  minimum 
Frobenius  norm  solution  to  AX  =  Y .  In  both  cases,  this  solution  is  A  =  YX+ .  Note  that  it  is 
possible  to  turn  an  under-determined  system  into  an  over-determined  system  by  truncating  the 
number  of  POD  modes  used  to  less  than  rri  (truncating  to  precisely  m  results  in  a  unique  solution 
when  the  data  is  full  column  rank,  with  no  loss  of  data).  A  least-squares  solution  of  this  form 
minimizes  the  error  in  Y ,  but  implicitly  assumes  that  there  is  no  error  in  X.  This  can  explain  why 
the  bias  in  DMD  (Eq.  (6.2.6))  is  dependent  on  Nx,  but  not  Ny  That  is,  in  the  least-squares  case 
DMD  can  be  viewed  as  finding 

A  :  Y  +  Ey  =  AX,  minimizing  ||£y||^, 

where  ||  •  \\y  denotes  the  Frobenius  norm  of  a  matrix.  When  doing  backwards-time  DMD  in 
Sect.  6.2.4,  we  conversely  assume  that  Y  is  known  exactly  and  minimize  the  error  in  X.  That 
is,  assuming  the  identified  dynamics  are  invertible,  we  find 

A  :  Y  =  A(X  +  Ex ),  minimizing  ||Ex||f- 

For  this  reason,  combining  forward-  and  backward-time  DMD  takes  into  account  the  error  in  both 
X  and  Y .  A  more  direct  means  of  doing  this  is  to  use  a  single  algorithm  that  finds  a  least-squares 
solution  for  the  error  in  both  X  and  Y .  It  is  possible  to  adapt  standard  TLS  algorithms  [53]  to  a 
DMD  setting,  which  we  perform  here.  We  seek 

r  Ex 

A  :  (Y  +  Ey)  =  A(X  +  E\),  minimizing  ||E||^,  where  E  = 

1  1  \py. 

The  expressions  Y  +  Ey  and  X  +  Ex  can  be  interpreted  as  Ym  —  Ny  and  Xm  —  Nx  ■  To  solve  for 
this,  we  can  rearrange  the  equation  to  obtain 

I4  -J][rt£v]=°'  (6'2'16) 
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We  would  now  like  to  assume  that  2 n  <  m.  This  might  not  be  the  case,  particularly  for  high¬ 
dimensional  fluids  data.  To  get  around  this,  and  improve  computational  tractability,  we  may 
project  Eq.  (6.2.16)  onto  a  POD  subspace  of  dimension  r  <  m/2,  to  obtain 


[A 


-I] 


X  +  Ex 
Y  +  Ey 


=  0. 


(6.2.17) 


This  POD  projection  step  is  in  contrast  to  the  TLS  DMD  formulation  in  [63],  where  a  projection 

~  X 

is  performed  onto  a  basis  determined  from  an  augmented  snapshot  matrix  Z  =  .  We  find  that 

the  present  formulation  yields  more  accurate  eigenvalues  in  a  number  of  examples.  Note  that  the 

P  +  can  have  rank 


nullspace  of  [ A  —  /]  is  r-dimensional,  meaning  that  the  2 r  by  m  matrix 
at  most  r. 


Y  +  Ey 


Let  the  full  SVD  of 


be  given  by  UEV* .  If  the  data  is  noisy,  we  should  expect  that  all  2 r 


diagonal  entries  of  E  are  nonzero.  By  the  Ec.kart- Young  theorem  [46],  the  nearest  (in  the  sense  of 
Frobenius  norm)  rank  r  matrix  will  be  given  by 


X  +  Ex 
Y  +  Ey 


=  UE1:rV 


where  Ei;r  contains  the  leading  r  singular  values  of  E,  with  the  rest  replaced  by  zeros.  We  then 
have  that 

X  +  EX] 

Y  +  Ey . 

where  Utj  are  r  by  r  matrices,  and  Vi  is  the  first  r  columns  of  V.  Rearranging  this  equation,  we 
obtain  the  total  least-squares  estimate  for  A: 


=  UZ1:rV*  = 


'Un 

U 1 2" 

Si  0 

~vf 

'Un^Vf 

U21 

U22_ 

0  0 

V* 
lv2  J 

U2iE1V1*_ 

A  =  U21UT1 


'11 


(6.2.18) 


Note  that  this  derivation  requires  that  U \  \  be  invertible.  While  the  derivation  includes  the  full  SVD 
of  the  augmented  data,  Eq.  (6.2.18)  indicates  that  we  only  need  the  first  r  columns  of  U ,  meaning 
that  only  a  reduced  SVD  is  required.  Algorithm  5  summarizes  this  total  least-squares  approach  to 
DMD,  which  we  refer  to  as  total  least-squares  DMD,  or  tlsDMD. 


Algorithm  5  (total  least-squares  DMD  (tlsDMD)). 

1.  Collect  data  X  and  Y,  and  project  onto  r  <  m/2  POD  modes  to  obtain  X  and  Y . 


2. 


Take  the  SVD  of 


X 

Y 


,  letting 


X 

Y 


U  ED. 


3.  Partition  the  2 r  by  2 r  matrix  JJ  into  r  by  r  sub-matrices,  letting  U 
only  the  first  r  columns  need  to  be  computed). 

4 ■  Compute  the  total  least-squares  DMD  matrix  A,  using  Eq.  (6.2.18). 


Uu 

U21 


U12 

U22 


(note  that 


5.  Compute  the  DMD  eigenvalues  and  modes  using  steps  4-5  of  Algorithm  2. 


An  alternative  and  more  focused  exposition  of  tlsDMD  is  given  in  [63].  We  note  that  Algorithm  5 
is  not  identical  to  that  presented  in  this  work  (due  to  the  lack  of  pre-truncation  of  POD  modes), 
however  we  find  that  Algorithm  5  gives  marginally  better  results  in  terms  of  the  accuracy  of 
identified  eigenvalues. 
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Figure  6.2:  (a)  Eigenvalues  (in  continuous-time)  identified  from  regular  DMD  (Algorithm  2,  dots) 
and  noise-corrected  DMD  (Algorithm  3,  crosses)  from  100  snapshots  of  data  from  Eq.  (6.3.1),  with 
At  =  0.1  and  crjj  =  0.01.  Only  one  of  the  complex  conjugate  pair  of  eigenvalues  is  shown.  The 
mean  and  95%  confidence  ellipse  of  1000  trials  are  given  for  each  data  set.  (b)  shows  the  mean  and 
95%  confidence  ellipse  for  the  same  data  set  for  Algorithms  2-5 


6.3  Results  with  synthetic  data 

In  this  section  we  will  test  our  proposed  modifications  to  DMD  on  a  number  of  examples.  Using 
known  dynamics  with  the  addition  of  random  noise  will  allow  us  to  examine  the  performance 
of  these  proposed  modifications  (Algorithms  3-5)  in  comparison  to  regular  DMD  (Algorithm  2). 
We  begin  by  considering  a  simple  2-dimensional  linear  system  in  Sect.  6.3.1.  In  Sect.  6.3.2,  we 
consider  the  same  system  with  an  expanded  set  of  observables,  which  tests  the  important  case  of 
high- dimensional  data  that  is  described  by  low-dimensional  dynamics.  Sect.  6.3.3  compares  the 
performance  of  the  proposed  modifications  of  DMD  to  other  DMD  variants  in  recent  literature, 
while  Sect.  6.3.4  considers  the  problem  of  identifying  dynamics  that  are  quickly  decaying  and 
obscured  by  dominant  modes  and  noise,  a  case  where  DMD-like  algorithms  could  be  of  most  use. 
Finally,  in  Sect.  6.3.5  we  analyze  how  the  proposed  DMD  modifications  treat  process  noise. 

6.3.1  Example:  A  periodic  linear  system 

We  consider  first  a  simple  two-dimensional  linear  system,  with  dynamics  given  by 

x  =  |  l  x'  (6.3.1) 

This  system  has  (continuous-time)  eigenvalues  Acq2  =  ±i,  so  gives  purely  periodic  dynamics, 
with  no  growth  or  decay.  We  discretize  with  a  timestep  At  =  0.1,  so  the  discrete-time  eigenvalues 
are  then  Aq2  =  e±Atl.  We  use  100  timesteps  of  data  (i.e.,  m  =  99),  corrupted  with  Gaussian 
white  noise  of  variance  cr2N  =  0.01.  The  identified  continuous-time  eigenvalues  from  both  regular 
DMD  (Algorithm  2),  and  the  direct  noise-correction  (Algorithm  3)  are  shown  in  Fig.  6.2(a),  for 
1000  different  trials  from  the  initial  condition  xo  =  [1  0.1] 2  .  We  assume  that  the  correction 
term  is  given  by  m.cr^In,  and  observe  that  this  corrects  almost  perfectly  for  the  bias  in  the  DMD 
algorithm  in  terms  of  identifying  eigenvalues.  Also  shown  in  Fig.  6.2(a)  are  ellipses  representing 
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Figure  6.3:  Error  in  the  estimated  propagation  matrix  A  arising  from  performing  DMD  and  ncDMD 
on  noise-corrupted  data  generated  from  Eq.  (6.3.1),  for  different  values  of  m  and  a2,.  In  (a)  the 
error  is  given  as  ||Mtrue  —  Mpred||F,  while  in  (b)  this  quantity  is  normalized  by  the  standard  deviation 
of  the  noise,  <7 at.  In  both  cases,  the  error  is  averaged  over  100  trials  for  each  m  and  a2N.  Note  that 
for  clarity,  (b)  excludes  the  two  largest  noise  levels  shown  in  (a) 


the  95%  confidence  region,  with  the  major  and  minor  axes  of  the  ellipse  aligned  with  the  principal 
component  directions  of  the  eigenvalue  data.  For  clarity,  in  the  presentation  of  subsequent  results, 
we  will  omit  individual  data  points  and  show  only  such  ellipses.  In  Fig.  6.2(b)  we  show  the  mean 
and  95%  confidence  ellipses  for  Algorithms  4  and  5.  As  with  ncDMD,  both  fbDMD  (Algorithm  4) 
and  tlsDMD  (Algorithm  5)  accurately  correct  for  the  bias  in  the  mean  of  the  identified  eigenvalue. 
Further  to  this,  fbDMD  and  tlsDMD  also  both  reduce  the  area  of  the  95%  confidence  ellipse,  which 
indicates  that  they  are  more  likely  to  attain  a  closer  approximation  to  the  correct  eigenvalue  on 
any  given  trial. 

Focusing  back  on  comparing  Algorithms  2  and  3,  we  show  results  for  a  variety  of  values  of  m 
and  <r)y  in  Fig.  6.3.  In  Fig.  6.3(a),  rather  than  looking  at  the  error  in  the  eigenvalues,  we  instead 
consider  the  Frobenius  norm  of  the  difference  between  the  true  and  identified  propagation  matrices, 
||Atrue-^pred||F-  For  very  small  noise,  the  correction  makes  little  difference,  since  the  random  error 
is  larger  than  the  bias  error.  For  larger  values  of  noise,  we  observe  that  the  error  saturates  when 
using  standard  DMD,  which  is  due  to  the  presence  of  the  bias  term  identified  in  Sect.  6.2.2,  which 
has  a  size  independent  of  the  number  of  samples,  m.  We  note  that  the  magnitude  of  this  bias  term 
is  proportional  to  a2N ,  as  predicted  by  Eq.  (6.2.10).  Evidence  of  this  error  saturation  phenomena 
can  also  be  seen  in  past  studies  of  the  effect  of  sensor  noise  on  DMD  [44,  164,  101].  After  this  bias 
term  is  corrected  for,  we  see  that  the  error  decays  proportional  to  m-1/2  for  all  values  of  noise, 
as  predicted  from  the  analysis  in  Appendix  1.  The  more  rapid  decay  in  error  with  m  for  small 
numbers  of  samples  seems  to  arise  from  the  fact  that  the  data  has  not  yet  completed  one  full  period 
of  oscillation.  Fig.  6.3  shows  the  corrections  to  DMD  made  using  both  the  sampled  (NxN^)  and 
theoretical  (mcr^I)  covariance  matrices.  Normally  the  sample  noise  covariance  would  not  be  known, 
and  so  we  demonstrate  here  that  the  theoretical  covariance  achieves  almost  the  same  decrease  in 
error.  Fig.  6.3(b)  shows  that  the  ncDMD  error  curves  collapse  when  the  error  is  normalized  by  the 
standard  deviation  of  noise,  ax  (note  that  we  could  also  multiply  the  error  by  the  SNR  to  get  the 
same  scaling). 

Fig.  6.4  shows  the  performance  of  Algorithms  4  and  5  on  the  same  data  as  Fig.  6.3.  Again,  we 
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Figure  6.4:  Error  in  the  estimated  propagation  matrix  A  arising  from  performing  DMD,  ncDMD, 
fbDMD,  and  tlsDMD  on  noise-corrupted  data  generated  from  Eq.  (6.3.1),  for  different  values  of  m 
and  a2N.  The  error  is  given  as  Untrue  —  ^4pred||.F>  and  is  averaged  over  100  trials  for  each  m  and  ct2n 


find  that  both  of  these  algorithms  can  prevent  the  error  saturation  present  in  standard  DMD,  and 
indeed  can  perform  noticeably  better  than  Algorithm  3  for  larger  noise  levels.  Algorithms  3-5  all 
appear  to  exhibit  the  same  asymptotic  behavior  as  the  number  of  snapshots,  m,  increases,  with  the 
error  decreasing  proportional  to  m~lt2 . 

A  common  means  to  mitigate  the  effect  of  noisy  data  is  to  collect  multiple  time-series  of  data, 
and  process  this  in  such  as  way  to  improve  the  results  over  just  using  one  data  set.  One  can  ask 
the  question  if  it  is  better  to  concatenate  the  snapshots  of  data  from  each  time  series  and  apply 
DMD  to  this  collection,  or  to  apply  DMD  to  phase-averaged  data.  Our  results  suggest  that  the 
latter  option  is  preferable  if  using  standard  DMD,  since  adding  additional  pairs  of  snapshots  will 
not  decrease  the  error  beyond  a  certain  level,  due  to  this  bias  saturation  at  large  m.  If  we  are 
using  ncDMD,  fbDMD,  or  tlsDMD,  however,  then  we  should  get  the  same  scalings  regardless  of 
which  option  is  chosen,  since  in  both  cases  the  error  should  be  proportional  to  p-1/2,  where  p  is 
the  number  of  trials  of  data  collected. 


6.3.2  A  periodic  linear  system  with  a  high-dimensional  state  of  observables 


This  example  considered  in  Sect.  6.3.1  has  m  n,  which  is  atypical  of  many  fluids  systems  for 
which  DMD  is  used.  To  consider  the  case  where  the  size  of  the  state  n  is  larger  than  the  number 
of  snapshots  m,  we  expand  the  state  of  our  system  to  include  time-shifts  of  the  data.  In  this  sense, 
we  have  new  observables  given  by 


Xfc 
Xfc- 1 


(6.3.2) 


X-k—q 


with  the  size  of  the  state  n  =  2(q  +  1).  This  periodic  system  can  equivalently  be  viewed  as  a 
traveling  wave,  which  is  now  observed  over  a  larger  spatial  domain.  Similar  data  (but  with  a  non¬ 
zero  growth  rate)  was  considered  in  [44]  and  [164],  Since  the  dynamics  are  still  only  two-dimensional 
despite  the  higher  dimensional  state,  we  use  only  the  first  two  POD  modes  of  the  data  to  identify 
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n  =  32 


n  =  256 


Re(Ac) 
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Figure  6.5:  Mean  and  95%  confidence  ellipses  of  continuous-time  eigenvalues  identified  by  applying 
regular  DMD  (Algorithm  2),  noise-corrected  DMD  (Algorithm  3),  forward-backward  DMD  (Al¬ 
gorithm  4)  and  total  least-squares  DMD  (Algorithm  5)  on  1000  trials  of  noisy  data  generated  by 
Eq.  (6.3.1),  and  observed  as  in  Eq.  (6.3.2).  Here  the  number  of  snapshots  rn  is  fixed  to  be  50, 
At  =  0.1,  and  =  0.1.  Only  one  of  the  complex  conjugate  pair  of  eigenvalues  is  shown 


a  2  x  2  propagation  matrix  A.  The  next  section  will  examine  alternative  means  of  performing  this 
dimensionality  reduction. 

Fig.  6.5  shows  the  statistical  results  (in  terms  of  DMD  eigenvalues)  of  performing  variants  of 
DMD  on  such  data,  using  m  =  50  and  a  range  of  snapshot  sizes,  n.  We  find  that  a  bias  exists  for 
regular  DMD,  but  the  magnitude  of  this  bias  decreases  as  the  size  of  each  snapshot  increases  (note 
that  the  scale  between  subplots  changes,  though  the  aspect  ratio  remains  the  same).  We  find  that 
Algorithms  3-5  all  outperform  regular  DMD  in  terms  of  giving  mean  (expected)  eigenvalues  that  are 
closer  to  the  true  value.  For  small  state  sizes,  Algorithms  4  and  5  also  also  give  a  smaller  confidence 
ellipse,  though  this  is  not  observed  for  larger  state  sizes.  As  the  size  of  the  state  increases,  the  bias 
component  of  the  error  of  DMD  (evidenced  by  the  difference  between  the  true  and  mean  identified 
eigenvalue)  becomes  smaller  relative  to  the  random  component  of  the  error  (indicated  by  the  size 
of  the  confidence  ellipse).  This  means  that  the  modifications  to  DMD  presented  in  Algorithms  3-5 
give  the  largest  improvement  when  the  size  of  the  state  is  small,  due  to  the  fact  that  in  this  regime 
the  bias  component  of  error  is  larger  than  the  random  component.  Note  that  these  conclusions 
may  be  predicted  from  the  scaling  laws  given  in  Eqs.  (6.5.2)  and  (6.5.5).  Moreover,  one  can  verify 
that  as  the  size  of  the  state  (n)  increases,  the  size  of  the  ellipses  decrease  proportional  to  vT1!2 . 

6.3.3  Comparison  to  other  modified  DMD  algorithms 

Without  any  modification,  applying  DMD  on  noisy  data  will  give  min(?n,  n)  eigenvalue-mode  pairs, 
many  of  which  may  be  mostly  or  entirely  due  to  noise,  particularly  if  the  underlying  dynamics  are 
low- dimensional.  For  this  reason,  a  number  of  modifications  of  DMD  that  aim  to  identify  a  small 
number  of  dynamically  important  modes  have  been  developed.  The  most  simple  means  of  reducing 
the  dimension  of  the  data  is  to  simply  project  onto  a  reduced  number  of  POD  modes,  which  is 
explicitly  mentioned  as  an  optional  step  in  Algorithm  2.  This  projection  step  was  also  used  within 
Algorithms  3-5  in  Sect.  6.3.2.  A  number  of  alternative  means  to  obtain  a  small  number  of  dynamic 
modes  from  DMD-type  algorithms  have  been  proposed,  as  briefly  mentioned  in  Sect.  6.1.  These 
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variants  all  start  with  the  observation  that  standard  DMD  can  be  formulated  within  an  optimization 
framework,  in  the  sense  that  it  identifies  a  least-squares  or  minimum-norm  propagation  matrix  for  a 
given  data  set.  [27]  proposes  a  modification  termed  optimized  DMD  that  seeks  to  find  optimal  low- 
rank  dynamics  that  best  matches  a  sequential  time-series  of  data.  While  the  fact  that  this  method 
optimizes  over  the  entire  time-sequence  of  data  rather  than  just  pairwise  snapshots  should  increase 
its  robustness  to  noise,  the  non-convexity  of  the  optimization  potentially  limits  its  utility.  Optimal 
Mode  Decompostion  (OMD,  [56,  164])  finds  an  optimal  low-dimensional  subspace  on  which  the 
identified  dynamics  reside,  rather  than  assuming  that  this  subspace  is  simply  the  most  energetic 
POD  modes.  This  approach  was  shown  to  give  an  improvement  on  the  DMD  eigenvalues  obtained 
for  noisy  data  in  [164],  Sparsity-promoting  DMD  (spDMD,  [71])  adds  an  l\  regularization  term 
that  penalizes  the  number  of  DMD  modes  with  non-zero  coefficients  in  the  approximation  of  the 
time-series  of  data. 

This  section  will  compare  OMD  and  spDMD  with  the  algorithms  presented  in  the  present 
work.  Of  the  algorithms  presented  here,  we  will  focus  on  fbDMD  (Algorithm  4),  which  was  found 
to  perform  equally  well  as  tlsDMD,  and  better  than  ncDMD,  in  Sects.  6.3.1  and  6.3.2.  Fig.  6.6  shows 
identified  eigenvalue  statistics  (mean  and  confidence  ellipses)  for  each  of  these  algorithms,  using  the 
same  data  as  that  for  Fig.  6.5.  We  observe  that  OMD  gives  a  more  accurate  mean  eigenvalue  that 
DMD,  and  a  confidence  ellipse  of  approximately  the  same  size.  spDMD  gives  a  mean  identified 
eigenvalue  that  is  closer  again  to  the  mean,  although  the  variance  in  the  eigenvalues  identified  for 
each  trial  is  larger.  We  note  that  spDMD  occasionally  produced  erroneous  results,  which  were 
excluded  as  outliers  from  the  statistical  analysis.  This  highlights  an  important  advantage  to  the 
modifications  to  DMD  presented  here  -  the  algorithms  are  given  in  closed  form,  and  do  not  rely  on  an 
appropriate  selection  of  parameters  and  tolerances  that  are  most  likely  required  for  an  optimization 
procedure.  In  all  of  the  cases,  fbDMD  (and  tlsDMD,  which  is  not  shown  but  barely  distinguishable 
from  fbDMD)  gives  the  best  estimate  of  the  true  eigenvalue. 

While  these  results  suggest  that  fbDMD/tlsDMD  is  more  accurate  than  OMD  and  spDMD,  we 
must  remember  that  the  results  from  one  data  set  do  not  show  the  global  superiority  of  any  given 
algorithm.  Indeed,  one  could  most  likely  find  data  sets  for  which  any  given  algorithm  is  superior 
(by  some  chosen  metric)  to  others.  We  conclude  this  section  by  noting  that  it  should  be  possible 
to  combine  the  optimization  procedures  presented  in  [27],  [164],  and  [71])  with  the  modifications 
to  DMD  presented  here.  Indeed,  a  simple  means  to  do  this  might  be  to  modify  Algorithm  4  so 
that  the  results  of  applying  a  given  algorithm  forwards  and  backwards  in  time  are  geometrically 
averaged,  as  in  Eq.  (6.2.15). 

6.3.4  Identifying  hidden  dynamics 

The  systems  considered  in  Sects.  6.3.1  and  6.3.2  could  be  considered  “easy”  in  the  sense  that  the 
dominant  dynamics  are  simple,  and  of  consistently  larger  magnitude  than  the  noise.  Indeed,  it  is  not 
difficult  to  qualitatively  identify  such  dynamics  by  eye  from  simply  looking  at  some  visualization  of 
the  data.  A  more  difficult  case  occurs  when  some  of  the  dynamics  are  of  low  magnitude  and/or  are 
quickly  decaying,  and  thus  might  quickly  be  lost  among  the  noise  in  the  measurements.  A  major 
benefit  of  data  processing  techniques  such  as  DMD  is  the  ability  to  identify  dynamics  that  might 
otherwise  remain  hidden.  With  this  in  mind,  we  now  consider  a  superposition  of  two  sinusoidal 
signals  that  are  traveling  across  a  spatial  domain  in  time,  with  the  amplitude  of  one  mode  growing, 
and  the  other  decaying: 

f(x,  t )  =  sin(/qa;  —  uqt)e7lf  +  sinfax  —  U2t)eY2t  +  na(x,  t ),  (6.3.3) 
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Figure  6.6:  Mean  and  95%  confidence  ellipses  of  continuous-time  eigenvalues  identified  by  applying 
regular  DMD  (Algorithm  2),  forward-backward  DMD  (Algorithm  5),  OMD  and  spDMD  for  noisy 
data  generated  from  1000  trials  of  data  generated  by  Eq.  (6.3.1),  and  observed  as  in  Eq.  (6.3.2). 
Here  the  number  of  snapshots  m  is  fixed  to  be  50,  At  =  0.1,  and  cr'2N  =  0.1.  Only  one  of  the  complex 
conjugate  pair  of  eigenvalues  is  shown 


where  we  set  k±  =  1,  c*q  =  1,  71  =  1,  &2  =  0.4,  uq  =  3.7  and  72  =  —0.2.  We  thus  have  the 
superposition  of  a  growing,  traveling  wave,  and  a  decaying  signal  that  is  quickly  hidden  by  the 
unstable  dynamics.  The  four  continuous-time  eigenvalues  of  this  data  are  71  ±  uq  and  72  ±  uq. 
This  data  is  again  similar  to  that  considered  in  [164]  and  [44],  if  we  neglect  the  decaying  dynamics. 
Fig.  6.7  shows  the  data  with  white  noise  of  standard  deviation  a  =  0.5.  Fig.  6.8  shows  the 
performance  of  various  DMD-type  algorithms  in  identifying  one  of  the  dominant  eigenvalues  (1  +  *) 
and  one  of  the  “hidden”  eigenvalues  (— 0.2  +  3. 7i).  Mean  eigenvalues  and  error  ellipses  are  computed 
from  1000  different  noise  samples.  Unsurprisingly,  all  methods  are  quite  accurate  at  identifying  the 
dominant  eigenvalue,  though  the  variants  proposed  in  the  present  work  show  improvements  in  both 
the  mean  and  scatter  over  the  1000  trials.  In  terms  of  the  hidden  eigenvalue,  we  observe  that  DMD 
(as  well  as  OMD)  estimates  a  decay  rate  that  is  almost  double  the  true  value.  In  contrast,  all  of 
ncDMD,  fbDMD,  and  tlsDMD  predict  the  eigenvalue  accurately,  with  a  reduction  in  the  error  of 
the  mean  eigenvalue  between  DMD  and  fbDMD  (for  example)  of  88%.  In  addition,  we  note  that 
the  scatter  in  the  identified  hidden  eigenvalue  across  the  trials  is  smaller  for  fbDMD  and  tlsDMD 
(as  indicated  by  smaller  confidence  ellipses). 

6.3.5  Differentiating  between  process  and  sensor  noise 

This  section  will  primarily  address  the  issue  of  comparing  and  distinguishing  between  the  effects 
of  process  and  sensor  noise.  We  consider  the  Stuart-Landau  equation,  which  has  been  used  as  a 
model  for  the  transient  and  periodic  dynamics  of  flow  past  a  cylinder  in  the  vortex  shedding  regime 
[95,  7].  In  discrete  time,  we  can  express  this  system  in  polar  coordinates  by 

rk+i  =rk  +  dt(nrk  -  rk  +  nr), 

9k+i  =  0k  +  dti'y  - /3rl  +  — ),  (6.3.4) 

rk 
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Position  0  0  Time 


Figure  6.7:  Visualization  of  data  generated  by  Eq.  (6.3.3),  with  k\  =  1,cji  =  1,  71  =1,  k-2  = 
0.4,  UJ2  =  3.7,  72  =  —0.2,  and  a  =  0.5 


Figure  6.8:  Mean  and  95%  confidence  ellipses  of  continuous-time  eigenvalues  identified  by  applying 
regular  DMD  (Algorithm  2),  OMD,  noise-corrected  DMD  (Algorithm  3),  forward-backward  DMD 
(Algorithm  4)  and  total  least-squares  DMD  (Algorithm  5)  to  1000  trials  of  noisy  data  generated 
by  Eq.  (6.3.3),  with  k\  =l,uj\  =  1,  71  =1,  &2  =  0.4, =  3.7,  72  =  —0.2,  and  a  =  0.5 
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Figure  6.9:  Eigenvalues  identified  using  (a)  DMD,  (b)  ncDMD,  and  (c)  tlsDMD  for  the  Stuart- 
Landau  equation  (Eq.  (6.3.4)),  with  100,000  snapshots  of  data  from  Eq.  (6.3.5),  with  ro  =  1, 
H  =  1,  7  =  1,  /3  =  0,  and  dt  =  0.01.  Data  with  sensor  noise,  process  noise,  neither  and  both 
are  considered,  with  noise  levels  for  process  and  sensor  noise  being  a2P  =  0.01  and  a2N  =  ICE4 
respectively.  Note  that  in  the  absence  of  sensor  noise,  DMD  and  ncDMD  are  identical 


where  we  have  included  process  noise  terms  nr  and  ng,  which  are  assumed  to  be  independent  in 
time,  and  sampled  from  separate  zero-mean  Gaussian  distributions  with  variance  dp.  We  take  as 
our  data  snapshots  of  the  form 


Xfc  =  ^e~Ji6k  e(-J+i)i0*.  ■■■  eJi9k]T , 


(6.3.5) 


for  some  integer  J .  We  may  add  sensor  noise  to  this  data  as  in  previous  sections.  For  fj,  >  0, 
Eq.  (6.3.4)  contains  a  stable  limit  cycle  at  r  =  yj Ji ,  with  period  27r/(7  —  /3/i).  Starting  on  the  limit 
cycle,  we  consider  data  with  process  noise,  sensor  noise,  neither,  and  both.  Without  any  noise, 
the  eigenvalues  identified  from  this  data  will  lie  upon  the  imaginary  axis,  at  locations  given  by 
Af,  =  ij(q  —  Process  noise  acts  to  perturb  the  system  from  its  limit  cycle,  which  ultimately 
leads  to  phase  diffusion,  and  a  “bending”  of  the  eigenvalues  such  that  they  instead  lie  on  a  parabola. 
The  behavior  of  this  system  with  process  noise  is  described  more  extensively  in  [8].  Fig.  6.9  shows 
the  results  of  applying  variants  of  DMD  on  data  generated  by  Eq.  (6.3.4)  with  fj  =  1,  7  =  1,  /3  =  0, 
and  dt  =  0.01,  with  data  collected  using  Eq.  (6.3.5)  with  J  =  10.  Applying  DMD  on  noise- free  data 
gives  eigenvalues  along  the  imaginary  axis,  while  data  from  the  system  with  process  noise  gives  a 
parabola  of  eigenvalues,  as  expected.  For  data  collected  using  Eq.  (6.3.5),  each  data  channel  will 
be  orthogonal  in  time,  and  will  contain  the  same  energy.  As  a  result,  sensor  noise  will  act  to  shift 
all  identified  eigenvalues  into  the  left  half  plane  by  the  same  amount,  as  observed  in  Fig.  6.9(a). 
Fig.  6.9(b)  shows  that  applying  ncDMD  accurately  corrects  for  this  shift,  for  the  system  with  and 
without  process  noise.  This  shows  that  it  is  possible  to  distinguish  between  the  effects  of  these  two 
forms  of  noise,  given  only  an  estimate  of  the  magnitude  of  the  sensor  noise.  That  is,  we  are  able  to 
eliminate  the  effects  of  the  noise  that  is  due  to  imperfections  in  our  observations,  while  retaining 
the  effects  of  actual  disturbances  to  the  system.  Fig.  6.9(c)  shows  that  tlsDMD  corrects  for  the 
effects  of  both  process  and  sensor  noise,  which  is  desirable  if  one  wishes  to  recover  the  dynamics 
of  the  noise-free  system.  The  results  for  fbDMD  are  not  shown,  but  were  very  similar  to  those  for 
tlsDMD.  The  ability  of  tlsDMD  and  fbDMD  on  process  noise  is  not  surprising,  since  they  treat  X 
and  Y  in  a  symmetric  manner,  and  thus  consider  phase  diffusion  both  forwards  and  backwards  in 
time. 
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6.4  Results  with  numerical  and  experimental  data 

Having  analyzed  the  performance  of  the  various  proposed  modifications  of  DMD  on  synthetic  data 
sets,  we  now  turn  our  attention  to  data  obtained  from  fluids  simulations  and  experiments.  We 
will  focus  on  the  canonical  case  of  the  unsteady  wake  of  a  circular  cylinder  exhibiting  periodic 
vortex  shedding.  In  Sect.  6.4.1  we  present  results  from  data  obtained  from  a  two-dimensional 
direct  numerical  simulation,  while  Sect.  6.4.2  considers  data  obtained  from  PIV  measurements  in 
a  water  channel. 

6.4.1  Cylinder  wake:  simulation  data 

We  use  an  immersed  boundary  projection  method  [139,  34],  with  a  domain  consisting  of  a  series 
of  nested  grids,  with  the  finest  grid  enclosing  the  body,  and  each  successive  grid  twice  as  large  as 
the  previous.  The  finest  grid  consists  of  uniformly  spaced  points  with  grid  spacing  equal  to  0.02D 
(where  D  is  the  cylinder  diameter),  extending  2D  upstream  and  4 D  downstream  of  the  center  of 
the  cylinder,  and  spanning  4 D  in  the  direction  normal  to  the  flow.  Each  successively  larger  grid 
contains  the  same  number  of  grid  points,  with  twice  the  grid  spacing  as  the  previous  grid.  The 
coarsest  grid  spans  24 D  in  the  streamwise  direction  and  16 D  in  the  normal  direction.  Uniform 
boundary  conditions  are  used  to  first  solve  the  Navier-Stokes  equations  on  the  largest  grid,  with 
each  smaller  grid  using  the  next  larger  grid  for  boundary  conditions.  The  numerical  scheme  uses 
a  3rd  order  Runge-Kutta  time  stepper,  with  a  time  step  of  0.01.D /Uoo,  where  U0 0  and  D  are  the 
freestream  velocity  and  cylinder  diameter,  respectively.  The  Reynolds  number  Re  =  U°^D  was  set 
to  be  100,  where  u  is  the  kinematic  velocity.  This  Reynolds  number  is  above  that  for  which  the  wake 
is  stable  (47,  [107]),  and  below  that  for  which  3-dimensional  instabilities  emerge  (approximately 
194,  [163]).  At  this  Reynolds  Number,  the  wake  is  hence  unstable,  and  approaches  a  single  periodic 
limit  cycle  characterized  by  a  von-Karman  vortex  street  in  the  wake.  The  data  to  be  analyzed  was 
taken  from  234  snapshots  of  the  vorticity  field,  spaced  0.1  D/U  time  units  apart.  This  corresponds 
to  approximately  4  complete  periods  of  vortex  shedding.  We  truncate  the  data  to  only  consider 
first  15  POD  modes.  These  first  15  POD  modes  contain  99.99%  of  the  total  energy  of  the  clean 
data,  and  92.96%  of  the  total  energy  of  the  data  after  the  addition  of  Gaussian  white  noise  with 
standard  deviation  a  =  0.2 U/D.  Thus  it  is  almost  entirely  noise  that  is  truncated  for  the  noisy 
data. 

Fig.  6.10  shows  results  from  applying  various  variants  of  DMD  to  such  data.  Though  not 
shown,  the  results  of  applying  tlsDMD  were  visually  indistinguishable  as  using  fbDMD.  Since  we 
are  artificially  adding  noise,  we  can  compare  the  results  using  noisy  data  to  those  generated  from 
the  noise-free  data.  When  applying  regular  DMD  to  noisy  data,  we  observe  significant  errors  in 
the  growth  rate  associated  with  the  highest-frequency  eigenvalues  (Fig.  6.10(a)).  For  an  oscillatory 
system  such  as  this,  the  DMD  eigenmodes  are  very  similar  to  the  POD  modes,  with  a  DMD  mode 
corresponding  to  Ac  ~  0  that  is  almost  the  mean  flow,  and  the  modes  associated  with  conjugate 
pairs  of  DMD  eigenvalues  corresponding  to  pairs  of  POD  modes  with  equal  energy,  see  [27]  for 
further  discussion  of  this  phenomenon.  This  means  that  the  observed  measured  eigenvalues  are 
in  line  with  the  analysis  given  in  Sect.  6.2.2  and  Appendix  1,  since  the  lower  energy  POD  modes 
oscillate  the  most.  We  can  see  the  effect  of  this  error  in  Fig.  6.10(b),  which  shows  the  prediction 
of  a  number  of  POD  coefficients  as  evolved  by  the  identified  system,  starting  from  the  true  initial 
condition.  The  dominant,  low  frequency  POD  modes  are  accurately  predicted,  but  the  higher 
“harmonics”  are  erroneously  predicted  to  decay  when  using  regular  DMD.  ncDMD  improves  the 
performance  marginally,  while  fbDMD  and  tlsDMD  both  almost  completely  remove  the  erroneous 
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Figure  6.10:  (a)  Eigenvalues  and  (b)  POD  coefficients  identified  from  applying  DMD,  ncDMD 
fbDMD,  and  tlsDMD  to  DNS  vorticity  data  from  a  cylinder  wake  at  Re  =  100.  Noisy  data  was 
corrupted  with  Gaussian  white  noise  with  a  =  0.2U/D 


decay  of  the  high-frequency  modes. 

As  well  as  considering  eigenvalues,  we  also  validate  in  Fig.  6.11  that  the  modifications  of  DMD 
do  not  adversely  affect  the  identified  DMD  modes.  This  is  shown  both  visually  in  Fig.  6.11(a),  and 
quantitatively  in  Fig.  6.11(b),  where  we  give  the  inner  product  (4>i, noisy,  4>i, clean)  of  the  ith  modes 
identified  from  clean  and  noisy  data,  where  we  have  pre-scaled  the  modes  to  be  of  unit  norm. 
We  enumerate  the  modes  by  the  imaginary  component  of  the  associated  eigenvalue,  with  mode 
0  corresponding  to  the  eigenvalue  on  the  real  axis.  For  modes  that  come  in  complex  conjugate 
pairs,  we  arbitrarily  consider  those  with  positive  imaginary  component.  We  see  that  both  fbDMD 
and  tlsDMD  marginally  outperform  regular  DMD,  in  terms  of  identifying  modes  that  are  at  least 
as  close  to  those  identified  from  noise-free  data.  The  decrease  in  the  inner  product  as  the  mode 
number  increases  is  indicative  of  noise  being  more  significant  in  higher-frequency  modes,  which 
contain  less  energy. 

6.4.2  Cylinder  wake:  experimental  data 

We  now  turn  our  attention  to  data  acquired  from  water  channel  experiment.  An  anodized  aluminum 
cylinder  of  diameter  D  =  9.5  mm  and  length  L  =  260  mm  was  immersed  in  a  recirculating,  free- 
surface  water  channel  with  freestream  velocity  Uoo  =  4.35 cm/s.  This  gives  a  Reynolds  number 
Re  =  DL^°°  =  413.  Further  details  of  the  experimental  setup  and  methodology  are  provided  in 
[148].  We  apply  variants  of  DMD  to  500  snapshots  from  a  vorticity  field  of  size  135  x  80.  Fig.  6.12 
shows  the  identified  eigenvalues  and  the  predicted  POD  coefficients  from  the  models  identified  from 
DMD  and  tlsDMD.  As  in  Sect.  6.4.1,  we  first  project  onto  the  15  most  energetic  POD  modes.  Note 
that  some  eigenvalues  (typically  quickly  decaying)  are  outside  the  range  of  the  plot.  As  was  the  case 
with  DNS  data,  we  observe  that  DMD  gives  eigenvalues  that  are  further  into  the  left  half  plane  that 
and  of  the  other  methods.  This  manifests  in  the  erroneous  prediction  of  decaying  POD  coefficients 
(Fig.  6.12(b)),  particularly  for  modes  that  are  less  energetic,  and  more  rapidly  oscillating.  We  thus 
conclude  that  more  accurate  low-dimensional  models  for  the  experimental  results  can  be  achieved 
by  using  tlsDMD.  We  note  that  this  improvement  can  be  attained  without  explicit  knowledge  of 
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Figure  6.11:  (a)  A  subset  of  the  DMD  modes  (real  components)  computed  from  applying  various 
variants  of  DMD  to  DNS  data  of  flow  around  a  cylinder,  (b)  Inner  product  between  (normalized) 
clean  modes,  and  modes  obtained  from  noisy  data  (with  cr/v  =  0.2Uoo/D) 


(b) 


POD  Mode  3 


Figure  6.12:  (a)  Eigenvalues  and  (b)  POD  coefficients  identified  from  applying  DMD  and  tlsDMD 
to  experimental  vorticity  data 


the  process  and  sensor  noise  characteristics. 


6.5  Discussion  and  conclusions 

It  was  shown  in  Sect.  6.2  that  simple  linear  algebraic  considerations  can  allow  us  to  derive  an 
estimate  for  the  bias  that  exists  in  all  standard  formulations  of  DMD.  This  subsequently  led  to 
the  formulation  of  the  three  modified  algorithms  that  we  suggest  can  be  used  to  eliminate  this 
bias.  Sect.  6.3  showed  that  this  predicted  bias  is  indeed  present  in  the  results  of  DMD.  Directly 
correcting  for  this  bias  term  (Algorithm  3,  ncDMD)  was  shown  to  almost  completely  eliminate 
this  bias.  While  this  modification  demonstrates  that  our  characterization  of  the  dominant  effects 
of  noise  was  accurate,  its  usefulness  is  limited  by  the  fact  that  it  requires  an  accurate  estimate  of 
the  noise  covariance.  Additionally,  the  presence  of  a  XT2  term  in  correction  factor  used  in  ncDMD 
makes  this  computation  unsuitable  for  cases  with  small  singular  values  that  are  not  truncated. 
On  the  other  hand,  the  correction  factor  in  Algorithm  3  may  be  applied  to  existing  DMD  results 
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with  minimal  computational  effort.  Algorithms  4  (fbDMD)  and  5  (tlsDMD),  which  do  not  require 
knowledge  of  the  noise  characteristics,  were  also  found  to  correct  for  the  bias,  and  also  were  able  to 
reduce  the  random  error  across  many  noise  realizations  (as  seen  by  smaller  associated  confidence 
ellipses  in  Fig.  6.2,  for  example).  Furthermore,  fbDMD  and  tlsDMD  were  found  in  Sect.  6.3.5  to  also 
compensate  for  the  effect  of  process  noise.  This  feature  could  be  desirable  or  undesirable,  depending 
on  the  purpose  for  which  DMD  is  being  applied.  Note  that  this  is  also  consistent  for  the  findings 
in  Sect.  6.4.2,  where  for  a  notionally  periodic  system,  tlsDMD  was  found  to  give  eigenvalues  very 
close  to  the  imaginary  axis,  despite  (presumably)  the  presence  of  both  sensor  and  process  noise. 

In  practice,  the  examples  examined  in  Sects.  6.3.4,  6.4.1  and  6.4.2  suggest  an  overarching 
principle:  while  regular  DMD  can  be  accurate  for  identifying  dominant  dynamics  that  have  much 
larger  amplitudes  than  the  noise  in  the  data,  accurate  identification  of  the  eigenvalues  associated 
with  lower  amplitude  modes  (and  in  particular,  their  real  components)  can  be  significantly  improved 
when  using  the  modified  DMD  algorithms  presented  here.  Conversely,  if  one  is  primarily  concerned 
with  the  identification  of  modes  and  their  frequencies  of  oscillation,  and  less  concerned  with  accurate 
identification  of  growth/decay  rates,  then  the  effect  of  sensor  noise  is  comparatively  minimal,  and 
subsequently  the  choice  of  DMD  algorithm  is  less  important. 

Fundamentally,  the  bias  in  DMD  arises  because  the  algorithm  is  essentially  a  least-squares 
algorithm,  which  is  designed  for  cases  where  the  “independent”  variable  (which  in  DMD  takes  the 
form  of  the  data  X )  is  known  accurately,  and  the  “dependent”  variable  Y  contains  the  noise/error. 
In  reality,  since  X  and  Y  should  both  be  affected  by  noise,  minimizing  the  error  in  both  the  X 
and  Y  “directions”  can  allow  for  a  more  accurate  answer  to  be  obtained.  One  drawback  of  tlsDMD 
is  that  it  requires  taking  the  SVD  of  a  larger  matrix.  Note  that  for  cases  where  n  >  m  (i.e., 
the  size  of  each  snapshot  is  larger  than  the  number  of  snapshots)  and  there  is  no  truncation  of 
POD  modes  corresponding  to  small  singular  values,  DMD  gives  the  minimum  Frobenius  norm  (of 
A)  solution  to  AX  =  Y.  In  this  case,  in  principle  neither  fbDMD  or  tlsDMD  should  yield  any 
improvements.  In  reality,  however,  if  there  is  noise  in  the  data,  then  we  do  not  necessarily  want  an 
exact  fit  to  the  data,  but  rather  an  unbiased  estimate  of  the  noise-free  dynamics.  We  may  obtain 
this  by  truncating  POD  modes  that  are  deemed  to  be  mostly  noise,  and  use  some  variant  of  DMD 
to  identify  the  remaining  dynamics.  tlsDMD  and  fbDMD  give  very  similar  results,  which  suggests 
that  fbDMD  can  be  viewed  as  a  computationally  cheaper  alternative  to  approximating  the  results 
of  tlsDMD.  Note  that  while  fbDMD  is  often  computationally  cheaper,  it  relies  on  being  able  to 
invert  the  matrix  Brn ,  which  might  be  an  ill-conditioned  operation  for  some  data. 

In  Sect.  6.3.3,  we  compared  the  variants  presented  here  with  two  recent  optimization  algorithms 
that  have  been  proposed.  The  results  show  our  algorithms  outperforming  both  sparsity  promoting 
DMD  and  OMD.  Note  that  since  these  algorithms  are  not  in  closed  form,  but  instead  contain  opti¬ 
mization  procedures,  the  results  depend  somewhat  on  the  specification  of  the  relevant  optimization 
parameters.  In  this  comparison,  our  use  of  Algorithms  3-5  relied  upon  the  projection  onto  a  low  di¬ 
mensional  subspace  before  applying  DMD-type  algorithms.  We  particularly  note  that  the  tlsDMD 
algorithm  proposed  here  is  slightly  different  from  that  given  in  [63]  due  to  this  POD  projection, 
which  we  found  empirically  to  give  improved  results.  We  suggest  that  this  is  because  the  initial 
truncation  of  low-energy  POD  modes  has  a  filtering  effect  that  better  isolates  the  true  dynamics, 
at  least  for  the  datasets  considered  here.  One  could  imagine,  however,  that  in  certain  cases  this 
projection  could  lead  to  significant  degradation  of  results.  For  example,  where  the  dynamically  im¬ 
portant  modes  are  highly  dissimilar  to  the  dominant  POD  modes,  the  flexibility  for  the  projection 
basis  to  be  modified  could  be  particularly  advantageous.  In  such  cases,  sparsity  promoting  DMD  or 
OMD  could  give  more  favorable  results.  In  general,  it  is  relatively  common  in  system  identification 
to  use  a  subspace  that  is  larger  than  the  dimension  of  the  underlying  dynamics,  and  then  later 
truncate  to  obtain  a  reduced  order  model  of  an  appropriate  size/rank.  This  can  be  particularly 
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important  when  the  dealing  with  specific  system  inputs  and  outputs  [112].  [73]  discuss  a  number 
of  ways  in  which  true  dynamic  modes  can  be  distinguished  for  noisy  modes,  in  the  context  of  the 
eigensystem  realization  algorithm.  [149]  further  discusses  how  DMD  modes  can  be  scaled,  from 
which  appropriate  modes  can  be  chosen.  The  spDMD  algorithm  in  [71]  essentially  automates  this 
procedure,  and  comes  with  the  additional  potential  advantage  of  not  requiring  a-priori  knowledge 
of  the  dimension  of  the  reduced  order  dynamics  to  be  identified.  Note  that  it  is  also  possible 
to  combine  the  modifications  to  DMD  proposed  here  with  the  OMD  and  spDMD  optimization 
procedures,  which  could  result  in  further  improvements  in  some  circumstances. 

Though  we  used  a  large  number  of  trials  when  testing  our  results  on  synthetic  data  in  order  to 
obtain  statistically  meaningful  findings,  in  reality  one  would  most  likely  not  have  this  luxury  with 
real  data.  In  this  case,  it  is  important  to  understand  for  the  size  and  quality  of  the  data  to  be 
analyzed,  both  the  best  algorithm  to  use,  and  the  amount  of  confidence  that  should  be  had  in  the 
results  of  the  chosen  algorithm. 

While  this  work  has  been  motivated  by  and  has  largely  focused  on  sensor  noise  (that  is,  noise 
which  only  affects  measurements,  and  not  the  system  dynamics),  the  characterization  and  removal 
of  process  noise  (i.e.,  disturbances  to  the  system  states)  is  entirely  another  matter.  Interestingly, 
the  effect  of  process  noise  was  identified  analytically  in  [8]  to  be  a  parabolic,  decay  in  the  growth  rate 
of  identified  eigenvalues  with  increasing  frequency.  It  turns  out  that  a  similar  effect  is  observed  here 
for  the  case  of  measurement  noise.  Isolating  sensor  noise  from  process  noise  (especially  with  limited 
prior  information  about  the  statistics  of  either)  is  an  important  and  challenging  task,  particularly 
when  dealing  with  more  complex,  turbulent  flows,  where  the  true  dynamics  exist  on  a  wide  range 
of  spatial  and  temporal  scales.  The  fact  that  DMD,  ncDMD  and  tlsDMD/fbDMD  each  perform 
differently  on  these  different  forms  of  noise  could  itself  be  an  important  tool  to  this  end. 

Particularly  in  experimental  data,  users  might  typically  preprocess  data  in  a  number  of  ways 
before  considering  applying  DMD-type  algorithms.  It  could  be  advantageous  to  investigate  precisely 
how  various  averaging  and  smoothing  operations  affect  the  subsequent  analysis  of  dynamics,  and 
subsequently  whether  such  post-processing  and  analysis  can  be  achieved  through  a  single  algorithm. 

Ultimately,  having  a  larger  selection  of  possible  algorithms  should  be  of  benefit  to  researchers 
who  desire  the  dynamical  information  that  DMD-type  algorithms  can  provide,  who  can  choose  based 
on  the  size  of  the  data,  amount  of  noise  present,  required  accuracy  of  the  results,  and  amount  of 
computational  resources  available.  One  of  the  major  advantages  of  DMD  (and  related  algorithms) 
advocated  in  [126]  is  the  fact  that  it  requires  only  direct  data  measurements,  without  needing 
knowledge  of  any  underlying  system  matrix,  thus  making  it  well  suited  to  use  on  experimentally 
acquired  data.  Inevitably,  however,  data  from  experiments  is  always  affected  to  some  extent  by 
noise.  It  is  thus  important  to  properly  understand  and  quantify  how  noise  can  influence  the  results 
of  DMD.  Conversely,  the  quest  for  high  quality  data  can  often  require  large  investments  of  both 
time  and  money.  Formulating  algorithms  that  are  more  robust  to  noisy  data  can  be  a  cheaper 
alternative  to  obtain  results  of  sufficient  accuracy.  As  it  becomes  easier  to  generate  and  store 
increasingly  large  datasets,  it  is  also  important  to  recognize  that  simply  feeding  larger  quantities 
of  data  (e.g.,  more  snapshots)  into  a  given  algorithm  does  not  guarantee  desired  improvements  in 
the  accuracy  of  their  outputs,  as  illustrated  in  Figures  6.3  and  6.4. 

The  problem  that  fluid  dynamicists  face  in  extracting  tractable  information  from  large  datasets 
is  not  unique  to  fluids,  and  rather  transcends  a  wide  variety  of  fields  of  study  (although  other  fields 
are  often  not  afforded  the  luxury  of  knowing  the  underlying  differential  equations).  It  is  valuable  to 
recognize  and  make  use  of  the  parallels  in  previous  and  current  developments  across  a  wide  range 
of  fields.  We  likewise  hope  that  other  areas  can  benefit  from  the  work  that  is  motivated  by  the 
desire  to  understand  how  fluids  flow. 
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Appendix  1:  Quantifying  the  size  of  the  bias  in  DMD 


We  seek  to  quantify  the  magnitude  of  this  bias  present  in  DMD  that  was  derived  in  Sect.  6.2.2,  sub¬ 
ject  to  certain  simplifying  assumptions  on  the  nature  of  the  data  and  noise.  If  the  noise  is  uniform, 
and  spatially  and  temporally  independent,  then  K(NXNX)  =  E (U*  NXNXU)  =  U*m,crxU  =  ma'frl, 
where  crx  is  the  variance  of  each  independent  component  of  the  noise  matrix.  Furthermore,  if  we 
assure  that  we  are  projecting  onto  the  POD  modes  of  the  clean  data,  then  (XX*)  =  E2,  where 
UEV*  is  the  singular  value  decomposition  of  X.  Thus  with  these  assumptions,  Eq.  (6.2.6)  can  be 
simplified  to  give 

E (Am)  =  A(I  -  maxY~2).  (6.5.1) 

The  (diagonal)  entries  E2  of  E2  have  the  interpretation  of  being  the  energy  content  of  the  ith  POD 
mode.  We  then  should  expect  that  E2  ~  mnqta2X:  where  cr\  is  the  RMS  value  of  the  elements  in 

E2 

the  data  matrix  X,  and  qt  =  Tracc^S2 )  is  the  proportion  of  the  total  energy  of  the  system  contained 
in  the  ith  POD  mode.  For  this  scaling,  we  make  the  assumption  that  adding/removing  rows  and 
columns  of  data  (i.e.,  varying  m  and  n )  does  not  affect  either  a\  or  The  bias  term  E-2  is 
a  diagonal  matrix  whose  ith  entry  has  a  size  ( e\,)i  proportional  to 


(e&);  ~ 


1 

nqiSNR 2’ 


(6.5.2) 


where  SNR  is  the  signal-to- noise  ratio.  Thus  sensor  noise  has  the  effect  of  reducing  the  diagonal 

~  2 
entries  of  the  computed  Am  matrix  by  a  multiplicative  factor  of  1  —  nqNai  ,  which  means  that 

POD  coefficients  are  predicted  to  decay  more  rapidly  than  they  actually  do.  This  effect  is  most 
pronounced  for  lower  energy  modes,  for  which  the  qt  is  smaller.  We  thus  expect  to  identify  with 
DMD  (continuous-time)  eigenvalues  that  are  further  into  the  left  half  plane  than  they  should  be  (or 
would  be  if  we  applied  DMD  to  noise- free  data).  [44]  argues  in  the  case  of  periodic  data  that  the 
growth  rate  of  the  eigenvalues  should  typically  be  the  most  challenging  to  identify,  since  there  are  a 
range  of  pre-existing  methods  that  can  identify  frequencies.  Here  we  have  argued  that  it  is  precisely 
this  growth  rate  that  is  most  affected  by  noise.  Importantly,  the  amount  of  bias  is  independent  of 
m ,  which  suggests  that  the  bias  component  of  the  error  will  be  particularly  dominant  when  we  have 
a  large  number  of  low- dimensional  snapshots.  Importantly,  this  suggests  that  one  cannot  always 
effectively  reduce  the  effect  of  noise  by  simply  using  more  snapshots  of  data,  since  the  bias  error 
will  eventually  become  the  dominant  error. 

While  we  can  now  quantify  the  magnitude  of  the  bias  in  DMD,  we  do  not  as  yet  know  how  it 
compares  to  the  random  component  of  the  error  that  would  arise  from  a  given  realization  of  noise. 
To  do  this,  we  will  estimate  the  typical  size  of  the  variance  of  individual  entries  of  A,  using  the 
standard  definition 


var 


(6.5.3) 


Referring  back  to  Eq.  (6.2.3),  if  we  exclude  terms  that  are  quadratic  or  higher  in  noise,  and  assume 
that  the  noise  covariance  matrix  is  sufficiently  close  to  its  expected  value,  we  find  that 


(YmX+)~  E 


(Y  +  Ny)(X  +  N x)(XX*  +  XN*x  +  NxX*  +  NxNx)-1  ~  YX+  -  E(1VX1V^)S 


YX+(XN*x  +  NxX*)  +  NyX*  +  YN*X 


E-2. 


Elements  of  the  terms  XNx,  NXY*,  NYX*,  and  YNX  are  uncorrelated  sums  over  m  random 
terms,  with  each  term  in  the  sum  having  variance  nqiaxax  where  as  before  qi  is  the  energy  fraction 
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in  the  ith  POD  mode.  This  means  that  the  sum  will  have  variance  rnnq1(j'\.(j\J .  Assuming  that 
YX+(=  A)  does  not  greatly  change  the  magnitude  of  quantities  that  it  multiplies,  and  assuming 
that  qi  remains  constant  when  varying  m  and  n,  this  means  that  we  find  that 


var 


a 


N 


mna 


x 


(6.5.4) 


Thus  the  expected  size  of  the  random  error  in  applying  DMD  to  noisy  data  is 

Cr  ~  m1/2n1/2SNR  <'6'5'5') 

Comparing  Eq.  (6.5.5)  with  Eq.  (6.5.2),  we  propose  that  the  bias  in  DMD  will  be  the  dominant 
source  of  error  whenever 

m1/2SNR  >  n1/2. 
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content  of  the  forces  analyzed.  In  addition,  the  flow  fields  are  studied  using  Dynamic  Mode  Decomposition 
(DMD). 

Next,  we  show  how  a  recently  developed  extension  to  DMD  can  be  used  to  obtain  nonlinear  reduced-order 
models  for  fluids  systems.  We  modify  the  extended  DMD  algorithm  to  include  a  Tikhonov  regularization 
step,  which  is  found  to  give  improved  results  for  the  purposes  of  nonlinear  system  identification.  The 
method  is  demonstrated  on  the  canonical  example  of  flow  past  a  circular  cylinder,  and  is  shown  to  be 
superior  to  classical  POD-Galerkin  projection 

Finally,  we  examine  the  effect  of  noise  in  DMD.  As  well  as  giving  an  explanation  for  a  previously  identified 
sensitivity  to  noisy  data,  three  variants  of  the  DMD  algorithm  are  presented,  all  of  which  perform  better  on 
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